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ABSTRACT 



This study explores the optimization method of simulated annealing for use in an 
analytic tool for counter-narcotics analysis. A model is developed, employing RAND 
Corporation’s Simulation of Adaptive Response model as an objective function evaluator, 
which optimizes interdiction asset locations relative to a sample smuggling network in the 
Caribbean region. In addition to asset location optimization, the response of the model 
to changing numbers of assets is also tested. Results indicate that this methodology has 
potential for use in the counter-narcotics program, and perhaps other network interdiction 
applications. Further research and testing are recommended. Military drug interdiction 
and the smuggling threat are discussed. 
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THESIS DISCLAIMER 



The reader is cautioned that computer programs developed in this research may not 
have been exercised for all cases of interest. While every effort has been made, within 
the time available, to ensure that the programs are free of computational and logic errors, 
they cannot be considered validated. Any application of these programs without 
additional verification is at the risk of the user. 
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I. INTRODUCTION 



A. BACKGROUND 

The Defense Authorization Act of 1989 designated the Department of Defense 

(DoD) "as the single lead agency of the Federal Government for the detection and 

monitoring of aerial and maritime transit of illegal drugs into the United States." The Act 

also tasks DoD with expanded support to civilian Law Enforcement Agencies (LEA’s) 

involved in drug interdiction. According to Secretary of Defense Dick Cheney, 

The Department of Defense is an enthusiastic participant in the nation’s drug 
control effort. We have significant resources at our disposal. We can make a 
substantial contribution to our national effort if we use our assets intelligently and 
efficiently. [Ref. 1: p. 4] 

To these ends, the President’s National Drug Control Strategy for 1991 estimated the 
DoD budget for drug interdiction at nearly $1.1 billion, more than double the actual 
expenditures for 1989. Obviously, the potential contributions by DoD are being taken 
very seriously. 

In 1988, the RAND Corporation developed the Simulation of Adaptive Response 
(SOAR) model to analyze the effect of military involvement in drug interdiction. SOAR 
is a stochastic network simulation of drug smuggling involving theoretical routes for drug 
shipments and the response of smugglers to being interdicted. Expanded involvement of 
the military in drug interdiction was evaluated by examining the effect that increased 
interdiction probability on the theoretical routes might have on the import price of illegal 
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drugs. The conclusions of the authors of the RAND study (Gordon Crawford, Peter 
Reuter, et al.) were not heartening. Although the limitations of the SOAR model caused 
it to be biased toward successful interdiction, it was found that "interdiction must be very 
stringent indeed to greatly affect U.S. drug consumption." [Ref. 2: p. 73] 
However, the model is simplistic in its approach to smuggling routes and the presence of 
interdiction forces. Routes are not associated with actual geographic areas due to a lack 
of available data, while interdiction assets are represented solely by a probability of 
interdiction on the generic routes. This simplicity leaves room for expansion of the model 
and further analysis. 

B. PROBLEM ANALYSIS 

Government-sponsored studies prior to that undertaken by RAND did not 
concentrate on the adaptation of smugglers, nor explicitly examine the potential for 
military contribution to the drug interdiction program [Ref. 3: p. 4], As part 
of a comprehensive study, SOAR was designed to aid in a policy analysis of the overall 
interdiction program. 

As a policy tool, SOAR did not require a geographic structure or specific 
representation of interdiction assets. For this reason, SOAR is not useful in operational 
planning or analysis. RAND’s conclusions, at a policy level, indicate that military 
involvement in drug interdiction may have limited effectiveness [Ref. 3: p. 130]. 
However, the Department of Defense has a major role in the program, and would benefit 
from the availability of an analytic tool which could be used to evaluate the effects of 
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operational planning and to explore, with a consistent measure, the effectiveness of 
different force and capability alternatives. 

Interdiction problems are frequently modeled as How networks for which 
optimization methods are used to determine which arcs or nodes should be cut or 
removed to provide an optimum result. In the case of a smuggling network, transit routes 
may be represented by arcs which are not well defined, or are completely unknown to the 
interdictors, making it difficult to entirely sever the arc. Even if it were possible, the 
smugglers themselves will eventually realize that the arc has been cut and will move to 
another route, adapting their methods to the interdiction. The SOAR model simulates this 
type of response on the part of drug smugglers, making it a fairly obvious choice to 
measure the effectiveness of interdiction on a smuggling network. 

Because routes used by smugglers may not be clearly defined, it is unrealistic to 
merely place assets on arcs representing segments of the routes. By developing a 
geographic structure for the routes, assets can be located by latitude and longitude, and 
the effects of assets on the routes (i.e., the probability of interdiction on routes) can be 
measured as a function of their range from the routes. This also would allow assets to 
affect multiple routes, for example, routes passing through a geographically restricted 
area, or chokepoint. Finding the optimum locations for assets relative to a smuggling 
network with stochastic flows will be difficult. As the locations of assets change, flows 
through the routes will also change, resulting in local and global optima which vary with 
the locations of assets. There are also likely to be many such optima, a result of different 
assets contributing the same probability of interdiction to routes or being located in 
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positions previously tested by other assets. In addition, the number of possible locations 
for assets is infinite, requiring some control over the locations tested if discrete 
optimization is to be used. This control must include geographically feasible areas, but 
should also reflect operational constraints, such as a requirement for certain assets to 
operate within a set range from each other or a base, and political constraints, such as 
foreign airspace or waters restrictions. 

An optimization model with the properties described in this analysis, which retains 
enough realism to be considered valid, would be of great benefit at the operational level 
of the military drug interdiction program. 

C. METHODOLOGY 

We will propose a solution to the need for an analytical tool which may be used for 
analysis and planning at the operational level of the military drug interdiction program. 
An exploration of the use of simulated annealing to optimize the location of counter- 
narcotics assets with respect to a drug smuggling network is presented. The objective 
function is an outcome of the modified SOAR model, in which the routes used by the 
smugglers are modeled with a geographical structure, and the physical location of 
interdiction assets (eg. ships, aircraft, aerostat balloons, observation posts, etc.) are used 
to determine the probability of interdiction in an area. From the SOAR model we extract 
the total cost to smugglers, a measure of effectiveness for one major goal of the overall 
drug interdiction program, reduction in consumption of illegal drugs. Interdiction raises 
the risks involved with drug smuggling, thereby increasing the costs associated with 
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personnel pay, equipment and drug replacement, etc. Eventually, this results in higher 
consumer prices for drugs. Economically, higher prices should cause a decrease in 
demand for drugs in the long run, that is, lower consumption by Americans. Increased 
crime among addicts to compensate for higher prices is not addressed. 

Use of SOAR as the objective function evaluator in a model which exhibits 
numerous local and global optimal asset locations eliminates the potential use of gradient- 
optimization methods. It is for this reason that we chose to apply simulated annealing 
as the optimization method to maximize the total cost to the smugglers. 

D. SCOPE 

This thesis briefly discusses the flow of cocaine from South and Central America 
into the southern United States, and the role the Department of Defense plays in 
interdiction efforts within this geographical area. 

The model presented is capable of evaluating the use of assets in the Atlantic, 
Pacific and Southwest Border regions; however, for the numerical results presented in this 
thesis, only routes and assets in the Caribbean and Gulf of Mexico are examined. Much 
of the information that is available about known smuggling routes and the assets involved 
in counter-narcotics operations is highly classified. In the interest of maintaining an 
unclassified example, some data used in the model is taken from the original SOAR 
model, while information on smuggling routes and counter-narcotics assets is intended to 
be representative only, and should not be construed as data representing actual 
circumstances. 



5 



II. THREAT AND INTERDICTION EFFORTS OVERVIEW 



A. CURRENT THREAT OF ILLEGAL DRUG TRAFFICKING 

A number of indicators show that progress has been made in the war against drugs. 
Overall drug use, drug-related medical emergencies, and student attitudes toward drug use 
have all declined past the goals set in the 1989 National Drug Control Strategy, shown 
in Table 1. [Ref. 4: p. 5] 



Table 1 NATIONAL DRUG CONTROL STRATEGY NATIONAL OBJECTIVES: Statistical 
significance ranges from 1 to 30%. 





GOAL 


ACTUAL 


CURRENT OVERALL DRUG USE 


- 10% 


- 11% 


CURRENT ADOLESCENT DRUG USE 


- 10% 


- 13% 


OCCASIONAL COCAINE USE 


- 10% 


- 29% 


FREQUENT COCAINE USE 


50% 

reduction in 
rate of increase 


- 23% 

actual decline 


CURREN T ADOLESCENT COCAINE USE 


- 20% 


- 49% 


DRUC.-RELATED MEDICAL EMERGENCIES 


- 10% 


- 18%' 


STUDENT ATTITUDES TOWARD DRUG USE 


- 10% 


- 28% 



Cocaine-specific indicators showed decreased purity and significantly higher prices 
in many areas of the country for the past year. Also, since 1985, "... reported current 
cocaine use (use in the past month) decreased some 72 percent." 
[Ref. 5: p. 1] 
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The effect on smugglers themselves may be deduced from the surrender of a 
number of Colombian drug lords and the level of effort smugglers put forth to get their 
product into the United States. Examples of this effort include the tunnel between 
Mexico and Arizona discovered in May 1990 [Ref. 6: p. Al] and the 
occasionally high tech hiding places used in vehicles, boats, and aircraft. 

On the other hand, many problems still exist which illuminate the threat to national 
security declared by President Reagan in 1986. Drug-related violence and murders remain 
at high levels within the United States, as does the amount of cocaine being produced 
throughout South America. Table 2 displays production levels for the three major 
producing countries from 1988 to 1990, the last year for which data was available. 



Table 2 COCAINE PRODUCTION ESTIMATE: Potential Cocaine HC1 Production by Country 
(tons). Sources: NNICC Reports, 1989 and 1990. 



COUNTRY 


YEAR 


1988 


1989 


1990 


PERU 


242 


440-468 


440-473 


BOLIVIA 


112 


264-458 


259-435 


COLOMBIA 


43 


77 


72 


TOTAL 


397 


781-996 


771-980 



Smugglers’ high tech methods display their increasing adaptability to interdiction 
efforts, as well as the effects of interdiction. Some smugglers are showing that they are 
willing (or crazy enough) to take on interdiction forces, as evidenced by incidents along 
the border with Mexico where smugglers have exchanged small arms fire with U.S. 
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military forces [Ref. 7: p. 1], and the attempt in April 1990 by smugglers 
linked to the Medellin cartel to buy Stinger missiles [Ref. 8: p. 3]. 

The evolution of narco- terrorism may be proof enough that we are indeed facing 
a war which must be fought by our armed forces outside the United States. Many 
revolutionary/terrorist organizations, such as M-19 and FARC in Colombia and Shining 
Path in Peru, are intimately linked to the drug trade, forming the basis for this narco- 
terrorism. These groups have shown decided anti-American leanings, and have used 
violence and threats against high level U.S. officials. One-time M-19 leader Ivan Mariano 
Ospina said in December 1984, "May these threats be carried out and may they be carried 
out in the entire world against the rapacious imperialism that lives at the cost of misery 
of exploited people...." [Ref. 9: pp. 36,37] Even state leaders have 

been involved, most recently illuminated by General Noriega in Panama, but also 
including Fidel Castro who has acknowledged the use of cocaine as a weapon, "We are 
going to make the people up there (the United States) white, white with cocaine [Ref. 9: 
p. 32]." 

Clearly, there are some signs of success in the drive to reduce illegal drug 
trafficking and abuse. However, it is also plain that much more effort is required in all 
phases of the war on drugs before the situation can be considered under control, let alone 
resolved. Involvement of the Department of Defense in reducing the flow of drugs into 
this country is only one part of the overall struggle, but, as Secretary Cheney said, the 
military can have a great impact in this effort. 
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B. WHAT IS DOD’S INVOLVEMENT NOW? 



Since the Department of Defense’s first, somewhat reluctant, involvement in 
counter-narcotics in the 1970’s, the military has been called upon for increasing support 
each year. The Defense Authorization Act of 1989 presented DoD with a primary role 
in the drug war. We now see all segments of the military accepting counter-narcotics 
assignments, including operational missions, such as ship and aircraft patrols with 
embarked law enforcement personnel, intelligence gathering and surveillance, and material 
support. The overall DoD chain of command for counter-narcotics is shown 
in Figure 1. 




Figure 1 DoD Counter-narcotics Chain of Command. 
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Obviously there is high-level attention being paid to the problem. The three Joint 
Task Forces command the military assets involved in actual day-to-day interdiction 
operations. 

1. Joint Task Force 4 

Located in Key West, Florida, JTF4 is in charge of operations in the Atlantic 
area of responsibility (AOR). As the first Joint Task Force, their mission included the 
creation of a joint fusion center for tactical intelligence and the communications necessary 
to assemble the data and provide a finished product to all users. JTF4 conducts their own 
detection and monitoring operations in the Atlantic AOR, and coordinates the operations 
of other agencies, such as the Coast Guard, Customs Service and local law enforcement 
activities. 

Although JTF4 has no dedicated assets, they do have tactical control of 
various personnel, ships, patrol and intercept aircraft, and intelligence collecting assets 
assigned from other DoD commands, the Coast Guard and Customs Service. 
[Ref. 10: p. 78] 

2. Joint Task Force 5 

Located in Alameda, California, JTF5 has a similar responsibility and mission 
in the Pacific AOR as that of JTF4 in the Atlantic. They conduct organic operations to 
detect and monitor aircraft and surface vessels suspected of smuggling, as well as 
coordinating other agencies’ operations within their AOR. JTF5’s mission includes 
integrating into the intelligence communications network and providing drug-related 
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intelligence to law enforcement agencies as appropriate. [Ref. 11] As is the case 
with JTF4, JTF5 has no dedicated assets, but takes tactical control of various assets 
assigned front DoD commands, the Coast Guard and Customs Service. 

3. Joint Task Force 6 

Located at Ft. Bliss, Texas, JTF6 has a much different mission from that of 
JTF’s 4 and 5. Specifically, "Joint Task Force Six plans and coordinates all DoD support 
requested by Federal, State and local law enforcement agencies within the Southwest 
Border region." [Ref. 12] 

Requests for support usually come through Operation Alliance, an agency 
comprised of representatives from Federal, State, and local law enforcement agencies of 
the Southwest Border states: Texas, New Mexico, Arizona, and California. Support 
provided includes observation posts, reconnaissance, terrain denial operations designed 
to deter smugglers from using areas along the border, training teams, military drug 
detection dogs, and construction. Personnel from all military services and various 
National Guard units have been called upon to perform these functions. 

C. WHERE MAY DOD’S INVOLVEMENT LEAD? 

There are as many different opinions concerning DoD’s future involvement in drug 
interdiction as there are people contemplating it. These opinions range from complete 
disassociation to complete involvement in all aspects, not just detection and monitoring. 
Many rational, and some not so rational, arguments are given for both extremes and 



everything in between. 



One of the main arguments against the use of military forces in counter-narcotics 
is the Posse Comitatus Act (Title 18 U.S. Code, Section 1385), which prohibits the direct 
use of military forces for civil law enforcement activities within the United States. 
Because of this law, military forces involved in drug interdiction operations may not 
perform search and seizure or make arrests. The military is limited to its monitoring and 
detection mission and support of law enforcement activities. 

Another strong argument is that drug interdiction operations have a detrimental 
effect on DoD’s primary mission to maintain a military capable of fighting a war to 
defend the U.S. or its interests. This argument contends that the fight against drugs is not 
truly a war, and that performing area patrols and ship boardings, detection and intercept 
of slow moving aircraft, and patrolling a border area are not suitable training operations 
for the military. Following the end of the Cold War (i.e., the reduction of the Soviet 
threat) drug smuggling takes its place among "principle threats ... that replace the old 
Soviet threat." [Ref. 13: p. 18] This means that drug interdiction 
operations need no longer be construed merely as training evolutions, but are actually part 
of the military’s primary mission. 

A third argument is the cost of military drug interdiction operations. It was 
estimated that the Air Force spent approximately $433,000 per drug bust in 1986. 
[Ref. 14: p. 7] Many believe that the $1 billion-plus military drug 
interdiction budget could be better spent on education and treatment programs. These 
programs are extremely important, since interdiction alone cannot halt the problem of 
drug abuse. However, education and treatment are more effective when drugs are more 
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expensive and less readily available. "Most leaders ol prevention and treatment programs 
recognize this; their task is made easier when drug enforcement works." 
[Ref. 15; p. 3] The goal of interdiction is to "complement and support our 

international drug control activities and domestic law enforcement programs, ... to create 
an integrated supply reduction program.” [Ref. 15; p. 65] Interdiction functions not only 
to seize drugs destined for the U.S., but to act as a deterrent to potential drug smugglers, 
and to force the price of illegal drugs to rise through increased risks the smugglers face. 
In addition, "Interdiction provides a highly visible sign to other nations of U.S. interest 
in reducing drug use." [Ref. 3: p. 2] This is a typical extension of the age-old Navy 
mission of showing the flag. Thus, the cost per seizure of military interdiction operations 
does not encompass the whole military contribution to the overall drug control strategy. 

Plainly, the military has an important role to play in combatting the drug problem 
faced by the United States. While no one wants to see the military take on civilian law 
enforcement roles, DoD has the manpower, assets, technology, and money to assist and 
supplement Federal, State, and local law enforcement agencies in reducing the supply of 
drugs. 

While DoD maintains its role in detection and monitoring of illegal drugs entering 
the country, there are many who are pushing for an even greater military involvement. 
In 1990, Senator Mitch McConnel (R-KY) proposed a plan for shooting down suspected 
airborne smugglers who refuse to land for inspection [Ref. 16: p. A19]. 
Although this plan was not put into effect, it did have a large amount of support. 
Military and LEA personnel are providing training and technical assistance for anti-drug 
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and counterterrorist forces in a number of production and transit countries. Operation 
Bahamas, the Turks and Caicos Islands (OPBAT), a multinational strike force in place 
since 1983, has been very successful in preventing smuggling flights into the Bahamas. 
It is certainly conceivable that similar strike forces could be established with other 
countries. 

The completion of the aerostat balloon network and the Caribbean Basin Radar 
Network, as well as proposed over-the-horizon radar systems, would significantly increase 
the ability to detect and monitor air traffic in smuggling areas, and act as a deterrent to 
smugglers. Aerostats are tethered radar-equipped balloons, flown at about 10,000 feet, 
which "provide surveillance of oceanic and land areas right down to the surface." 
[Ref. 17: p. 7], with a range of approximately 200 nautical miles. The 

Caribbean Basin Radar Network is a series of ground based radars designed to provide 
nearly complete coverage of the Caribbean basin. Figure 2 displays the proposed 
coverage of the aerostat network [Ref. 18: p. 84J. 

RADM Walter Leland, USCG, said, "We know that when the aerostats are flying 
the traffic stops." [Ref. 19: p. 7] However, there are widely differing 
reports on the operational reliability of the aerostat balloons. Some Air Force and 
Customs Service officials call the system too expensive and fragile, claiming the balloons 
are not operational more than 60% of the time because of severe weather and 
maintenance problems [Ref. 20: p. 15]. At the same time, Pentagon 
figures showed the balloons to be operational 70-75% of the time. 
[Ref. 21: p. SI 4] Such discrepancies cast doubt on the effectiveness of 
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Figure 2 Aerostat Radar Network. 



these systems, in turn making it harder to get funding to continue and expand the 
program. 

Finally, there have been a number of popular novels recently published which 
provide more speculation about the possibilities of the military’s future role in drug 
interdiction. Most notable of these are Tom Clancy’s Clear ami Present Danger and Dale 
Brown’s Hammerheads. 

Although many widely ranging proposals for military involvement in counter- 
narcotics and numerous arguments against any use of the armed forces in drug 
interdiction exist, with the current trend toward downsizing the armed forces it is difficult, 
if not impossible, to predict exactly what DoD’s future role in the nation’s drug control 
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strategy will he. It is fairly certain however, that the U.S. military will continue to be 

involved. A recent report by Rep. Les Aspin (D-Wis.), Chairman of the House Armed 

Services Committee, foresees the continuation of a military role: 

The drug problem must be attacked from both sides -- both the demand for drugs 
and the supply of drugs must be reduced.... Although the United States will 
continue to rely primarily on nonmilitary means in its effort to curtail the 
international trade in drugs, Americans will also want U.S. military forces engaged 
in the struggle with drug traffickers, primarily in an interdiction role. 

[Ref. 13: p. 17J 



16 



III. SIMULATION MODEL 



A. THE ORIGINAL SOAR MODEL. 

1. Background. 

The Simulation of Adaptive Response (SOAR) model is part of a RAND 
Corporation study on the drug interdiction program, undertaken at the request of the 
Under Secretary of Defense for Policy in February 1986. This study culminated with the 
publication of analysis results in three RAND reports. 

The primary conclusion of RAND’s effort was that: 

Increasing drug interdiction efforts are not likely to greatly affect the availability 
of cocaine in the United States.. ..This conclusion is driven primarily by the small 
share of total drug distribution costs that are accounted for by the smuggling sector. 
Only about 10 percent of the final price of cocaine comes from smuggling costs and 
profits. [Ref. 3: p. xi] 

This conclusion does not take into account other reasons for military involvement in drug 
interdiction, such as those discussed in the previous chapter. It does, however, suggest 
that the need exists for a tool to be used in counter-narcotics analysis. RAND’s SOAR 
model provides a reasonable starting point for developing such a tool. 

2. Model Description. 

SOAR is a dynamic network simulation which allows drug smugglers a 
number of air, sea, and land smuggling routes. These routes are generic in that they are 
not associated with particular geographic routes. The route used for a particular shipment 
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of drugs is probabilistically selected, depending on the costs of shipping drugs via that 
route and the expected risk costs for shipping via the method associated with that route 
based on a perceived probability of the shipment being seized. This probability is 
estimated by a time-weighted history of unsuccessful shipments on that route, and will 
change with time as shipments succeed or are interdicted. 

Time between shipments is an exponential random variable with an input 
mean. The quantity of each shipment is an input constant, as are the capacity and costs 
of the routes and shipping methods. SOAR does not attempt to model the adaptation of 
changing shipment sizes in response to interdiction severity. The capacity of the route 
and the size of the shipment are used to determine the number of trips that each shipment 
will require. These trips occur on the same route at the same time, leading to a saturation 
factor which increases the probability that the shipment is seized. 

Interdiction in SOAR is modelled by a general probability of interdiction for 
each route. "Military resources are treated. ..simply as means for augmenting particular 
activities, thereby raising the probability of success in those activities." [Ref. 3: p. 8] 
This probability is an input value and can be changed for predetermined time periods 
(phases). When multiple trips are required, the probability of interdiction on the route, 
for that shipment, is increased as a result of the saturation factor. 

The smugglers’ overall goal is to ship as large a quantity of drugs as possible 
at their perceived lowest cost. Costs to the smuggler include the cost of the drug, 
operating costs for the method of shipment, the cost of replacing lost assets if a shipment 
is seized, and pay to personnel. Smuggling personnel are assumed to be risk averse, 
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therefore pay is assumed to change as the square of the perceived risk involved in using 
a particular route. Increasing interdiction rates will affect the overall cost to the 
smugglers by raising each of these costs. 

For their model, the RAND authors attempted to estimate: 

• Quantities shipped, by route, in a given year. 

• Number of shipments, by route, in a given year. 

• Number of vessels identified as suspicious, by route, in a given year. Of those, the 
number pursued, by route in a given year. Of those, the number resulting in 
seizures, by route, in a given year. 

• Estimates of the compensation resulting from the likelihood of prison. 

• Estimates of smugglers’ nonrisk compensation and profits. (Ref. 2: p. 9] 

As can be expected, little concrete data was available, resulting in the use of informed 
guesses for many of the input parameters. 

Output of the model includes an echo of the input parameters and the statistics 
described in Appendix A, Sections A.2. and B. 1 . respectively. The output is an average 
of an input number of runs desired. Ten runs were used in the RAND analysis, providing 
"an adequately precise estimate of the overall means, especially compared with the 
imprecision of some of the input data and some of the assumptions incorporated in the 
model." [Ref. 2: p. 24] This output was not used in our analysis, but is included in the 
appendix for information. 

Computational details of the SOAR model are described in Chapter IV of 
Simulation of Adaptive Response, and are reasonably easy to interpret from the model 
itself. 
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B. THE SOAR SIMULATION AS THE OBJECTIVE FUNCTION EVALUATOR. 



1. Measures of Effectiveness. 

In evaluating drug interdiction efforts, many measures of effectiveness have 
been proposed. Civilian law enforcement agencies are concerned with numbers of 
seizures and smugglers prosecuted. News reports regularly include the street value of 
drugs seized. Because of its unique role of detection and monitoring and support to law 
enforcement, and the fact that the armed forces are prohibited from direct participation 
in the arrest and seizure phases of interdiction, the military’s measure of effectiveness 
(MOE) is based on the number of supporting missions provided. This MOE is not very 
useful for a model which attempts to find optimal locations for interdiction assets. 
Although it indicates the level of military support to reducing the How of drugs into the 
United States, it provides no quantitative measure of the effectiveness of the overall drug 
interdiction program. 

The quantity of drugs seized is a popular MOE, as it provides an evaluation 
of how much drugs are prevented from entering the country. However, it is possible that, 
when interdiction efforts are most successful, no drugs at all would be seized. That is, 
shipments could be deterred, smugglers could be forced to use the most expensive 
methods of shipping to avoid interdiction, etc. 

The major goal of drug interdiction is to reduce the consumption of illegal 
drugs. Military assets are an integral part of increasing the effectiveness of interdiction, 
thereby helping to achieve this goal. In addition to seizing drugs and equipment, 
interdiction increases the risks which smugglers face and causes them to change their 
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methods of operation. All of these effects of interdiction result in higher costs to the 
smugglers, costs which are eventually passed on to the consumers of illegal drugs. 
Although some users are likely to be insensitive to increased drug prices, "the body of 
economic support for a relationship between price and consumption is too strong to allow 
much doubt that, at least in the long run, higher prices will lead to lower consumption." 
[Ref. 3: p. 20J Therefore, if it can be measured or estimated with some confidence, cost 
to the smugglers is an MOE which reflects the overall goal of the drug interdiction 
program. 

As computed by the SOAR model, total cost to the smugglers includes the 
costs of unsuccessful shipments, thereby incorporating the MOE of the amount of drugs 
seized. Parameters for an alternate objective function value could easily be extracted 
from the simulation, or the optimization algorithm could be applied to any desired 
objective function. 

The simulation used as the objective function evaluator for the optimization 
model presented in this thesis is a slightly modified version of the original SOAR model. 
The changes involve the attempt to model smuggling routes with a geographical structure 
and calculating the probability of interdiction based on the location of military assets. 

2. Model Assumptions. 

In any model, many simplifying assumptions are required. This is especially 
true in attempting to model a process of this scale. Perhaps the biggest assumption lies 
in the use of the SOAR model outcome as the objective function to be optimized. That 
is, we assume that the simulation provides a reasonable estimate of the total cost to the 
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smugglers. This includes accepting as reasonable the assumptions made in constructing 
SOAR, which are described in the RAND publications. 

Another major assumption involves modelling the drug smuggling routes. 
Except for the natural chokepoints in the Caribbean, little is known about the exact routes 
being used. Therefore, the routes presented in this model are assumed to be 
representative of area routes of unknown width. In the Mona Passage, Yucatan Channel 
and other chokepoints, the width of sea routes is constrained by the presence of land on 
either side. However, in the Pacific, Southwest Border area, and for many air routes, the 
width could be hundreds of miles. Assets which provide coverage of a route are assumed 
to be patrolling the area which that route represents. Figure 3 displays a general idea of 
smuggling routes as interpreted by the Drug Enforcement Agency. 




Figure 3 Drug flows from South America into the United States. 
Source: DEA Map 1 - 1987 (Revised 1989) 
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A third assumption is made when computing the probability of detection on 
the routes. Each asset represents a sensor which adds a probability of interdiction to the 
route that it covers. It is also assumed that an asset covering a leg of the route acts 
independently of other assets on that and other legs, and that detection on any leg of a 
route can be taken as an event independent of all other legs on that route. Realistically, 
a patrol aircraft on a route, which detects a suspicious surface contact is likely to be in 
either direct or indirect contact with a surface asset on that leg or a subsequent leg, giving 
information about the presence of that contact to the surface asset and increasing the 
surface asset’s probability of detecting and interdicting that contact. Similarly, a surface 
asset may detect a suspicious aircraft and report it, etc. For modeling purposes, however, 
we believe these assumptions to be reasonable. Independence between assets and routes 
will result in SOAR providing an upper bound on the measure of effectiveness. This 
consequence should be considered when estimating the probabilities of interdiction input 
to the model for each asset type against each method of smuggling. 

3. Data Structures and Computations. 
a. Smuggling Routes. 

The representation of routes is accomplished using a hierarchical 
adjacency list of geographical points. The geographic points are listed in a data File 
similar to that shown in Appendix A Section A.5., which includes points of origin for 
drug shipments, intermediate route points, and points of entry into the United States. 
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The information for the routes themselves is taken from another data file. 



A description of the inputs and an example of the data file are in Appendix A Section 
A. 3. 

b. Assets. 

Appendix A Section A.4. describes the input information and a data file 
for interdiction assets. Although we have listed assets as SHIP1, PAIR1, etc., asset need 
not necessarily refer to a single ship, aircraft or ground unit. An asset may be considered 
whatever mix of elements is required to achieve the range and probability of interdiction 
input to the model. For simplicity, we will continue to refer to assets as single units. 
Asset types are used to describe assets with the same range and probability of interdiction. 
By assigning additional asset types, the matrices containing this information could be 
expanded to include unique input data for particular assets. 

c. Computation of route probabilities. 

A shipment of drugs will be interdicted on a particular route with a 
probability computed in the subroutine PROBCOMP. In this subroutine, the legs of each 
route are checked to determine if assets are located such that the leg falls within the range 
of the assets. Up to this range, the asset contributes a constant probability of interdiction 
to the route. Beyond this range, the probability decays exponentially as a function of the 
additional distance. The probability of interdiction on a route is then computed as simply 
the combination of independent probabilities for all assets covering all legs of the route. 
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The probability decay function chosen results in probabilities for 
individual assets decreasing to approximately 13.5% of the original value at twice the 
asset range. For example, an asset with a given constant probability of 0.10 to 100 miles 
will have a probability of interdiction equal to 0.0135 at 200 miles. This method of 
computing probabilities can be supported by considering the given range of an asset as 
that area which can be patrolled with a constant probability, and exponential decay of the 
probability representing the decreasing effectiveness of the patrolling as the distance 
between the asset and the route increases past that range. 
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IV. THE ANNEALING MODEL 



A. GENERAL DESCRIPTION 

The simulated annealing method is so-called for its resemblance to the annealing 
of solids, in which a material is heated to the point where all particles arrange themselves 
randomly in the liquid phase, and is slowly cooled so that all particles arrange themselves 
in the low energy ground state of a corresponding lattice. [Ref. 22: p. 2] 

The method starts at either a predetermined or random initial point. In the feasible 
space of the independent variables, a random walk samples the objective function 
[Ref. 23: p. 210] (analogous to the energy state of a physical system). 
Each step providing an improved objective function value is automatically accepted. 
Steps which are detrimental to the objective function value are accepted according to a 
Boltzmann probability function, p = exp(-AE/K B T), which is dependent on the energy and 
temperature of the physical system. This allows the path to walk out of local and global 
optimal points [Ref. 23: p. 211]. The process is often called a Metropolis loop, named 
for the author who used a model simulating the evolution of a solid to thermal 
equilibrium in 1953 [Ref. 24: p. 14]. Obviously, the algorithm could run 
indefinitely, walking in and out of the global optima. One way to prevent this from 
happening is to specify the number of steps the algorithm is to take. However, this in no 
way guarantees reaching the global optimum. Recently much research has been put into 
developing cooling or annealing schedules. Essentially, such schedules decrease the 
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probability of accepting a detrimental step by decreasing the temperature ol the system 
as the global optimum point is approached. The random walk is terminated when a 
predetermined number of successive trials does not produce an acceptable objective 
function value. [Ref. 23: p. 212] 

B. APPLICATIONS 

The simulated annealing method has been applied to optimization problems in 
numerous fields, including physics, operations research, numerical analysis, biology, 
materials science, game theory, code design, etc. [Ref. 22: p. 14]. 

There are many advantages and disadvantages associated with the simulated 
annealing method. The algorithm is applicable to many different optimization problems 
and is easy to implement. The greatest disadvantage is the potentially prohibitive amount 
of time required to converge to a near-optimal solution. This is dependent on the cooling 
schedule and the step size used in the random walk. [Ref. 25: p. 91] 

C. THE GENERALIZED ANNEALING ALGORITHM 

As stated previously, simulated annealing algorithms have been used almost 
exclusively in applications which minimize the objective function. The probability 
function is especially well-suited to cases where the objective function minimizes to zero, 
or where the optimal value is known. In these cases the probability of accepting a non- 
improving step is driven to zero as the optimum is approached. 

In the case where the optimal value is not known, the typical algorithm calls for 
selecting some initial estimate of the optimal value (<|) m ) and allow the random walk to 
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sample the objective function value (<f>), continuing until <J> - 4> In becomes negative, then 
decrease <j) m repeating this process as required. On the other hand, if the initial estimate 
of <)) m is too low, it is increased as part of the cooling schedule, so that cp - m will 
approach zero. [Ref. 23: p. 213] 

We have attempted one of many possible applications of this method to a 
maximization problem, by taking an initial guess for the optimum as an input value, as 
well as the number of steps to be taken before reducing the guess. Cooling the guess 
(<p m ) is accomplished by keeping track of the highest value found so far (<p B ), and 
reducing <J) m according to the function <(> m+1 = <f> m - oc * ((p m - <p B ), where a is a 
predetermined fraction which controls the cooling speed. This ensures that <|> ni is never 
reduced past the best value found in the random walk. If the random walk produces a 
value which is greater than either the current <)>„, or the input guess, then both are 
increased a significant amount above the higher value (we have chosen <]) m = 1.5 * <p B ). 
The cooling of <f> m again proceeds according to the previously described method. 

In our application of the algorithm, we begin with a large step size (Ar) in order to 
sample a significant portion of the feasible region. As <J> m is reduced (or cooled) we also 
reduce the step size to a minimum value (Ar min ) by a similar equation, sampling a tighter 
area of the region. When a potential point of convergence is found, Ar is reduced to a 
small value (Ar low ) in an attempt to find an acceptable move in the very close vicinity. 
If an acceptable move is found, the step size is reset to the previous value (Ar old ), and the 
annealing continues. 
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This algorithm is not guaranteed to locate the objective function value associated 
with a global optimum. The value (3 in the probability function controls the rejection 
ratio, and will also determine, in part, how close to the optimal value the algorithm 
converges. For very low values of (3, the algorithm will converge at, or very close to, the 
optimum, when the optimum value is known. 

The feasible space and step generation used for our model are very simple, 
consisting of upper and lower limits on the latitude and longitude for asset locations, and 
random steps with size Ar. For example, we did not attempt to prevent assets from 
moving onto land, restricted airspace, or within search range of other assets. These are 
the more detailed geographic, operational, and political constraints, such as those 
discussed in Chapter I, Section B, which could be used to control possible asset moves. 
Such constraints could be represented in a subroutine that tests feasible locations and a 
set of rules controlling the possible moves of assets relative to each other. 

Our algorithm, shown in Figure 4, is very similar to that presented by Bohachevsky, 
Johnson and Stein for a minimization problem allocating ballistic missile interceptors 
[Ref. 26]. The following terms used in the algorithm are defined: 

• <\> is the objective function value. A subscript of 0 indicates the current value, 1 
indicates the value at the new location. Other subscripts are as defined above. 

• Q is the feasible space, that is, the set of latitudes and longitudes where assets may 
be located. 

• x is the set of latitudes and longitudes for all mobile assets. A subscript of 0 
indicates the current location, 1 indicates the new location. 

• Ar is the step size. Subscripts min and low indicate the minimum step size used 
when cooling and the step size desired in attempt to find a more precise solution. 
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Step 1 : Compute the objective function value for initial locations, <(> 0 . 

Set 4> m = input guess. 

Step 2: For each mobile asset, generate two independent standard normal random variates, 
A, and A 2 , and compute the components of the random direction vector d: 
d, = A/(A,**2 + A 2 **2)*‘0.5, i = 1,2. 

Step 3: Compute candidate locations for assets, set x, = Xo + Ar * d. 

Step 4: If the candidate locations are infeasible, ie, x g Q, repeat Step 2. 

Step 5: Compute the new objective function value, <)>,, and set A<}> = <t>, - <|> 0 . 

Step 6: If the step increases the objective function value (A<J) > 0), accept: 
set k = 0, Xq = x,, <(> 0 = (p , , and n = n + 1 . 

If the step size has been reduced in Step 8 (k > S), reset step size Ar = Ar old . 

If <j>i > <j> 0 , set <)> e = <]),. 

6a: Increase <j> m and initial guess if new objective function value is greater, ie, 

lf <t>) > <t>m. Set 4>m = 15 * 4>1- 

If <J) , > initial guess, set guess = 1.5 * <)>,. 

Set n = 0. 

6b: If N accepted steps taken (n > N), cool <)> m and the step size: 
set <L = 4>m * « * (<L ' <t>e) and Ar = Ar - a * (Ar - Ar m J. 

Set n = 0. 

Go to Step 2. 

Step 7: If the step is a decrease (A<J> < 0), compute probability of acceptance and 
determine whether to accept or reject. 

Set p = exp( p * Atp /(<}> m -<t> 0 )). Generate a uniform random variate V. 

If V > p, reject, set k = k + 1 . If V < p, accept, set k = 0, Xq = x 1( <t> 0 = <J> t . 

Step 8: If the stopping rule has not been reached (k < S), go to Step 2, 

otherwise, if k = S, attempt to find a more precise solution with small step size, 
save current step size, Ar old = Ar, set Ar = Ar low , and go to Step 2. 

Step 9: If the final stopping rule has not been reached (k < 2 * S), go to Step 2, 
otherwise stop and return <J> 0 as the optimal solution found. 



Figure 4 Annealing algorithm to optimize interdiction asset locations relative to a 
stochastic network. 
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• P is the control parameter for the probability function (analogous to Boltzmann’s 
constant). 

• a is the fraction for reducing <{) m and Ar 

• N is the number of accepted steps at each value for <|) n) 

• S is the stopping rule, the control for the number of attempted moves rejected 
before the solution is declared found. 

• k is used to count the number of consecutive rejected moves 

• n is used to count the number of accepted moves for each level of <f» m 
The FORTRAN code for the model can be found in Appendix B. 

D. RESULTS 

The optimization model developed in this thesis was applied to a smuggling 
network of eleven air and sea routes in the Caribbean region. Figure 5 displays this 
network, formulated from the sample data sets in Appendix A. Some areas of the routes 
used are intuitive, for example, those running through the Yucatan Channel and the Mona 
and Windward Passages. However, these are meant to be illustrative of actual routes, and 
are not based on intelligence or other official information. The model was run on a 
number of different computers, SUN SPARC 1+ and 2 workstations, Tektronix 4336, 
Silicon Graphics IRIX, and an AMDAHL 5990-500 mainframe computer. The four 
UNIX systems yielded identical results, while the VMS-based mainframe computer found 
similar results to the UNIX systems. 

Two applications of the model were tested. The first examined the locations at 
which the annealing algorithm converged. Ten runs each were performed with five 
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mobile assets, three ships and two aircraft. The resulting locations from one run arc 
plotted in Figure 6. It is clear that some of the locations could easily be improved. This 
is due in part to the value of the probability function control parameter, p, used to ensure 
convergence of the algorithm in a reasonable amount of time. Each of the ten runs found 
different final locations, but all appeared reasonable and the optimal values found were 
similar, with mean equal to 75153.2 (dollars x 1000) and a standard deviation of 7204.1. 

The second application explored the response of the model to a changing number 
of assets. Each run used equal numbers of ships and aircraft. Parameters of assets are 
from the sample data set in Appendix A. We found that the function. 
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<J> = 615430.6 - 71.07 * (80 - ASSETS) ** 2.07, determined using non-linear regression, 
fit the data well (R 2 = 0.9922). This particular function is dependent on the number of 
assets at the last data point, and is extremely simplistic, with total numbers of assets as 
the independent variable. Runs with more than 80 assets will change the parameters of 
this function, but will not significantly increase the maximum. Although the fit of the 
curve is good, the variance at each level of assets cannot be assumed equal, mainly due 
to the small sample size at the levels (2). 

The regression model itself is not a fundamental result, but can be used to obtain 
useful information. With 80 assets we are very nearly saturating the routes, obtaining 
little additional benefit from adding assets. We will assume for the following discussion 
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that 80 assets induce the highest possible total cost to the smugglers. In Figure 7 we have 
plotted a curve depicting the fraction of maximum value achieved versus assets. This 
fraction represents the regression function value divided by 615430.6, the function 
evaluated at 80 assets. From this, we see that 50% of the maximum attainable effect on 
the smugglers is obtained with 24 assets, 75% with 36 assets. In other words, 75% of the 
maximum effectiveness is achieved with fewer than half the maximum number of assets. 
It is highly unlikely that such quantities of assets would be available for employment at 
any particular time. However, these results are entirely dependent on the specific 
parameters used and our simple asset mix. 




Figure 7 Fraction of maximum value found (at 80 assets) vs assets 
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V. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

In this thesis we have demonstrated that simulated annealing can be employed to 
find near-optimal locations for interdiction assets relative to a network with stochastic 
flows that are determined by a simulation. Our approach also demonstrated an application 
of this type of model to explore the relationship between the capabilities of interdiction 
assets, their number, and the resulting effect on the overall effectiveness of the 
interdiction effort. Although our model uses a simple representation of interdiction assets, 
combinations of assets, and the geographic, operational and political constraints which 
would normally be used to restrict feasible asset locations in the drug interdiction 
problem, we believe that the results demonstrate this methodology’s viability for 
development into the type of analytic tool described in Chapter I, Section B. 

The primary obstacle to further use of this methodology in counter-narcotics 
operational analysis and planning is the questionable reliability and limited amount of 
information available on smuggling activities. The credibility of results from this type 
of model, including the validity of SOAR, is highly dependent on the accuracy of 
information used in constructing the route network, and assigning costs to smugglers and 
parameters to assets. The suspect reliability of smuggling data will make it difficult to 
assess this credibility, short of actually employing operational assets and observing the 
resulting effects on smuggling and drug use. Even this level of verification would not 
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necessarily provide an accurate measure of credibility since the targets of drug interdiction 
are the only true source of the information required, and they are not likely to come 
forward with such details. 

B. AREAS FOR FURTHER RESEARCH 

The next logical step in the development of an analytic tool useful to the counter- 
narcotics program is to apply the methodology explored in this thesis, or the model 
presented itself, to an expanded problem. Such an expansion should include all primary 
cocaine smuggling regions (Atlantic/Caribbean, Pacific, and Southwest Border), 
constraints reflecting the geographic and political considerations for asset locations, and 
operational inputs, such as budget feasible asset combinations and capability restrictions 
of assets. These constraints could easily be included in a subroutine which tests proposed 
asset moves for feasibility. Analysis of the results from an expanded problem will 
provide a more realistic evaluation of the utility that may be derived from this type of 
tool. 

Expanding the size of the problem produces a major drawback in the use of the 
model developed in this thesis, extensive computer run times. With a network of only 
eleven routes, the five asset runs described in our results averaged seven hours each on 
the SUN SPARC 1+ workstation. Adding routes to the problem will significantly 
increase the run time of SOAR. Furthermore, our testing used only one trial of SOAR 
to obtain each new objective function value, while future testing should use the average 
of multiple trials to reduce the variability of the simulation results. This will increase run 
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times even more, by approximately a multiple of the number of trials. For these reasons, 
analysis of the SOAR model to optimize performance, or translation into a more efficient 
programming language, and the use of high speed computers would be beneficial to future 
efforts. 

A final consideration involves the variety of measures of effectiveness (MOEs) 
being used by drug interdiction agencies, as discussed in Chapter III, Section B.l. 
Development of a well-defined, quantifiable MOE for the entire drug interdiction program 
would be very valuable to future analysis. The MOE produced by the SOAR model for 
our optimization, total cost to the smugglers, relies heavily on information which may not 
be particularly accurate, as would be the case for any MOE that reflects the goals of the 
interdiction program as a whole. Continued intelligence efforts to gather accurate 
information on illegal drug smuggling activities and improvements in the timeliness of 
such data will help in defining an overall MOE and assessing the credibility of results 
from analytic tools used in counter-narcotics planning, and will make the drug interdiction 
program more effective. 
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A. MODEL INPUT AND OUTPUT 



A. INPUT DATA SETS 
1. ANNEAL DATA 
a. Variables 



Variable name 


Description 


STEPSIZE 


Step size for the annealing 
algorithm in miles. Real. 


STEPSIZEMIN 


Minimum step size for the 
algorithm in miles. Real. 


STEPSIZELOW 


Step size in fine search 
for optimum value in miles. 
Real. 


BETA 


Control parameter for 
annealing. Adjusted to 
control the 

acceptance/rejection ratio. 
Real. 


S 


The stopping rule for the 
algorithm. This number of 
consecutive rejections stops 
the algorithm. Integer. 


ALPHA 


Parameter used to decrement 
the estimate of the optimal 
value and step size. 

Adjusted to control the speed 
of cooling. Typically between 
0.01 and 0.2. Real. 
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GUESS 



Inital estimate of the 
optimal objective function 
value. Real. 



B1GN 



b. Sample data set 



The number of steps to be 
taken at each level of 
cooling. Integer. 



INPUT PARAMETERS FOR THE ANNEALING ALGORITHM 



STEPSIZE (MILES) 
MIN STEP SIZE 
SMALL STEP SIZE 
CONTROL PARAMETER 
STOPPING RULE 
ALPHA 

GUESS AT MAXIMUM 
STEPS PER LEVEL 



75.0 

40.0 
5.0 

1.5 
100 
0 . 1 

100000000.00 
1000 



2. SOAR DATA 
a. Variables 



Variable name 



Description 



NUMMETHODS 



Number of methods for 
smuggling drugs. Integer. 



NUMPHASES 



Number of phases. Integer. 



NUMROUTES 



Number of routes available. 
Integer. 



NUMTRIALS 



ENDTIME 



Number of times the analysis 
period is to be simulated. 
Reported results will be 
averaged over the trials. 
Integer. 

Number of days to be 
analyzed, not including the 
runin period. Integer. 
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RUNIN 



IX 



Number of days to run the 
simulation before commencing 
the analysis. Integer. 

Seed for random number 
generator. Integer. 



DRUGCOST 



The cost at the source of a 
kilogram of cocaine. Real. 



EXSHIPMENTINTVL Expected time between 

shipments. Real. 

EXSHIPMENTSIZE Expected shipment size. Real. 



METHODNAME(M) 



The name of drug smuggling 
method M. Character. 



CAPCOST(M) 



CAPACITY(M) 



The cost to the smuggler of 
an unsuccessful shipment 
using method M. Does not 
include the cost of the drug 
or the cost associated with 
the route. Real. 

Maximum amount of drug that 
may be shipped by method M. 
Real. 



b. Sample data set 



METHODS 

3 

DRUGNAME 

COCAINE 

METHOD 

AIR 

SEA 

LAND 

01 

61 



PHASES ROUTES TRIALS RUNIN ENDTIME 
2 3 2 120 60 

DRUGCOST SHIPMEAN SHIPSIZE DAILY AMT 



7500.0 0.71 

RISK COMP RISK EXP 

1200000.0 2.0 

1600000.0 2.0 

10000.0 2.0 



250.0 350.0 

SEIZE COST MAX SHIP 
200000.0 2000.0 

40000.0 16000.0 

5000.0 50.0 



MEMORY 
0 . 1 



SEED 

7243 
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3. ROUTES DATA 



a. Variables. 



Variable name 


Description 


ROUTENAME(R) 


The name of route R. 
Character. 


ROUTECOST(R) 


The cost to the smuggler for 
using route R, incurred 
whether or not the shipment 
is successful. Real. 


ROUTEMETHOD(R) 


The index of the smuggling 
method that is used on 
route R. Integer. 
(l=Air,2=Sea,3=Land). 


ROUTELOC(R) 


The index of the geographic 
area of route R. Integer. 

( 1 =Atlantic,2=Pacific, 
3=Southwest Border). 


FIRST(R) 


The index of the first point 
of route R. Integer. 


NEXT(R,P) 


The index of the point 
following P on route R. 
Integer. 
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EN 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 



b. Sample data set . 



ROUTENAME 


ROUTECOST 


ROUTEMETHOD 


POINTS 










MONA PASS 
1 10 


16000 .00 
11 38 




2 


1 


WINDWARD 
2 12 


16000 . 00 
13 14 


37 


2 


1 


YUCFLA1 
3 15 


16000 . 00 
36 




2 


1 


YUCALA1 
3 15 


16000.00 

35 




2 


1 


YUCTEX1 
3 15 


16000.00 

32 




2 


1 


YUCFLA2 
4 15 


16000.00 

36 




2 


1 


YUCFLA3 
4 15 


16000.00 

35 




2 


1 


YUCTEX2 
4 15 


16000.00 

32 




2 


1 


CUBA 
3 36 


20000.00 




1 


1 


YUCALA2 
7 34 


20000.00 




1 


1 


YUCTEX3 
7 31 


20000.00 




1 


1 



ROUTELOC 
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4. ASSET DATA 



a. Variables. 



Variable name 


Description 


NUMASSETS 


Total number of assets, 
mobile and fixed. Integer. 


NUMASSETTYPES 


Total number of asset types. 
Integer. 


ASSETNAME(A) 


The name of asset A. 
Character. 


LATASSET(A) 


The latitude of asset A. 
Real. 


LONASSET(A) 


The longitude of asset A. 
Real. 


ASSETTYPE(A) 


Index of the type of asset A, 
for example, 

1 =air,2=ship,3=land, 
4=aerostat balloon. Integer. 


ASSETLOC(A) 


Index of the geographic 
location of asset A. Same as 
ROUTELOC. Integer. 


ASSETRANGE(AT.M) 


Range for asset type AT 
versus smuggling method M for 
which the probability of 
interdiction is a constant. 

Real. 


IPD(AT,M) 


Probability of interdiction 
for asset type AT versus 
smuggling method M on a 
route, which is constant out 
to range ASSETRANGE. Real. 
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b. Sample data set. 



TOTAL NUMBER OF ASSETS = 8 

NUMBER OF ASSET TYPES = 4, PATROL AIR , SH I PS , LAND UNITS, AEROSTATS 



ASSET 


NAME 


LATITUDE 


LONGITUDE 


TYPE 


LOCALE 


1 


PAIR1 


16.00 


75.85 


1 


1 


2 


PAIR2 


17.45 


80.95 


1 


1 


3 


SHI PI 


19 . 00 


68.00 


2 


1 


4 


SHIP2 


19 . 00 


75.00 


2 


1 


5 


SHIP3 


21 .98 


86.18 


2 


1 


6 


CUDLOE KEY EL 


24.41 


81.30 


4 


1 


7 


MARFA TX 


30.19 


104.01 


4 


3 


8 


GREAT EXUMA 


23.33 


75.47 


4 


1 



ASSET RANGE 


VS METHOD 








METHOD 




ASSET 


AIR 


SEA 


LAND 


AIR 


100.0 


100.0 


25.0 


SHIPS 


50.0 


50.0 


0.0 


LAND 


10 . 0 


0.0 


5.0 


BALLOONS 


2 00.0 


0.0 


0.0 


P{ INTERDICT) 


' ASSET VS 


METHOD 






METHOD 




ASSET 


AIR 


SEA 


LAND 


AIR 


0.1 


0 . 1 


0 . 1 


SHIPS 


0.1 


0.1 


0.0 


LAND 


0.1 


0.0 


0.1 


BALLOONS 


0.1 


0.0 


0.0 



5. GEOGRAPHIC DATA 



a. Variables. 



Variable name 



Description 



POINTNAME(P) 



The name of geographic 
point P. Character. 



LAT(P) 



The latitude of point P. 
Real. 



LON(P) 



The longitude of point P. 
Real. 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 
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b . Sample data set . 



POINTNAME LATITUDE LONGITUDE 



PT GALL I NAS 


12 


.15 


71 , 


.45 


RIO HACHA 


11 . 


.30 


73 , 


.00 


SANTA MARTA 


11 , 


.15 


74 , 


.10 


TURBO 


8, 


.00 


76 , 


.45 


BUENAVENTURA 


3 , 


.58 


77 , 


. 50 


TUMACO 


1 , 


.45 


78 , 


.45 


MEDELLIN 


6 , 


. 10 


75, 


.47 


BOGOTA 


4 , 


.32 


74 , 


.15 


CALI 


3 . 


.20 


76, 


.33 


SMONAPASS 


17 , 


.37 


68 , 


.00 


NMONAPASS 


20 , 


.50 


68. 


. 00 


SWINDWARD 


17 , 


.95 


75 , 


.18 


NWINDWARD 


20 , 


.50 


73 , 


.82 


STRAIGHTS 


22 


.72 


76 , 


.91 


YUCATAN 


21 , 


.98 


86, 


.18 


CENT MEXICO 


25, 


. 00 


103. 


.00 


NOGOLES AZ 


31 , 


.20 


110. 


.55 


BISBEE AZ 


31 , 


.30 


109. 


.55 


DOUGLAS AZ 


31 , 


.20 


109. 


.30 


YSLETA TX 


31 , 


.42 


106. 


. 18 


PRESIDIO TX 


29 , 


.33 


104 , 


.23 


DELRIO TX 


29 , 


.21 


100. 


.52 


DOLORES TX 


27 , 


.42 


99 , 


.47 


BROWNSVILLE TX 


25, 


.95 


97 . 


.30 


SAN DIEGO CA 


32 , 


.41 


116 , 


.57 


CALEXICO CA 


32 , 


.41 


115. 


.30 


??1 NM 


31 . 


.40 


108 , 


.50 


??2 NM 


31 , 


.40 


108 , 


.25 


FREEPORT TX 


28, 


.56 


95. 


.21 


GALVESTON TX 


29 , 


.18 


94 . 


.48 


HOUSTON TX 


29 , 


. 46 


95, 


.21 


PORT LAVACA TX 


28, 


.36 


96 , 


.38 


SABINE TX 


29 , 


.44 


93 , 


.54 


MOBILE AL 


30 , 


.42 


88, 


. 03 


PENSACOLA FL 


30 , 


.25 


87 , 


.13 


ST PETE FL 


27 , 


.47 


82 , 


.38 


FLAMINGO FL 


25, 


. 10 


80 , 


.55 


DAYTONA BCH FL 


29 , 


.11 


81 , 


.02 


MIAMI FL 


25 , 


.45 


80. 


.11 


LOS ANGELES CA 


34 , 


.00 


118 . 


.15 


BAJA TIP 


22 , 


.64 


109 , 


.82 


S MEXICO 


15, 


.53 


96. 


. 18 


PACIFIC1 


13 , 


.32 


105 , 


.36 


PACIFIC2 


30, 


.00 


118 , 


.00 


PACIFIC3 


20 , 


.40 


107 , 


.27 


PACIFIC4 


15. 


.50 


105 , 


.00 
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B. OUTPUT DATA SETS 



1. SOAR OUTPUT 

Output of the model includes an echo of the input parameters, and the 

following statistics: 

1. Expected attempts per trial. The average number of shipments attempted during the 
period being analyzed. 

2. Expected successes per trial. The average number of successful shipments during the 
analysis period. 

3. Expected interdictions per trial. The average number of unsuccessful shipments 
during the analysis period. Expected attempts = Expected successes + Expected 
interdictions. 

4. Success rate. The proportion of shipments that were successful. 

5. Interdiction rate. The proportion of shipments that were unsuccessful. 

6. Cost of incomplete shipments (in thousands). The average cost to the smuggler 
because of unsuccessful shipments, including the cost of the method (such as an 
airplane), the cost of the drug, and the cost of the route (such as gasoline). 

7. Cost of completed shipments (in thousands). The average cost to the smuggler 
because of successful shipments, including the cost of the drug and the cost of the 
route. 

8. Total cost to smugglers (in thousands). The sum of the cost of incomplete shipments 
and the cost of completed shipments. 

9. For each drug, the quantity that the smugglers attempted to ship. 

10. For each drug, the quantity that the smugglers successfully shipped. 

11. For each drug, the quantity that the smugglers lost because of unsuccessful shipments. 

12. For each route, the expected attempts, successes, and failures are reported. 

13. For each phase and each route, the expected attempts, successes and failures are 
reported. 
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14. For each drug and each route, the expected quantities shipped, captured, etc.|Ref. 2: 
pp. 42,43] 



2. ANNEALING MODEL OUTPUT 



MAX TOTAL COST TO SMUGGLERS = 78047.95 

(IN THOUSANDS) 

OPT VALUE FOUND IN 6630 STEPS 

831 STEPS NOT ACCEPTED 

4195 INFEASIBLE MOVES FOUND 



FINAL LOCATIONS 

ASSET LATITUDE LONGITUDE 



PAIR1 


22 


.715 


84 , 


.834 


PAIR2 


23 , 


.539 


77 , 


.752 


SHI PI 


20 


.728 


83 , 


.596 


SHIP2 


20 


.333 


86, 


.613 


SHIP3 


21 


. 687 


68 , 


.579 


CUDLOE KEY FL 


24 


.410 


81 , 


.300 


MARFA TX 


30 


. 190 


104 , 


.010 


GREAT EXUMA 


23 , 


.330 


75. 


.470 



FINAL PROBABILITY OF INTERDICTION BY ROUTE 



ROUTE 


ROUTENAME 


PROBINTERDICT 


1 


MONA PASS 


0 . 1774 


2 


WINDWARD 


0.2130 


3 


YUCFLA1 


0.3206 


4 


YUCALA1 


0.3176 


5 


YUCTEX1 


0.3176 


6 


YUCFLA2 


0.2910 


7 


YUCFLA3 


0.2879 


8 


YUCTEX2 


0.2879 


9 


CUBA 


0.2249 


10 


YUCALA2 


0 .2909 


11 


YUCTEX3 


0 .1692 


CPU TIME 


= 17704, 


.3 SECONDS 


MAX VALUE 


: FOUND = 


7 . 83 9 17 E+07 


FINAL GUESS VALUE = 


8 . 10845E+07 


NEXT LOC SEED = 298087429 

RUN COMPLETED 
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APPENDIX B. FORTRAN CODE 



PROGRAM ANNEAL 

WRITTEN BY LT JAMES J. HENRY IV, USN 
NAVAL POSTGRADUATE SCHOOL, CODE 360 
OPERATIONS ANALYSIS CURRICULUM 



PROGRAM VARIABLES 

INCLUDE ' COMMON . FOR ' 

REAL PROBMOVE, 

Sc TLAT (MAXASSETS) , TLON (MAXASSETS ) , 

Sc OB J V AL 0 , OB J VAL 1 , OB J V ALM , OB J VALB , DELOB J , 

L V, P,TOTCOST, CPUTIME, 

Sc A ( 2 ) , DENOM , TARRAY { 2 ) 

INTEGER K, LITTLEN, I , ROUTENR , MOVES , NOTMOVE, INFEAS 
LOGICAL FEASIBLE, ACCEPT, FINISH , GOODINPUT 

CALL SUBROUTINE TO READ IN THE INITAL DATA 



CALL GETDATA (GOODINPUT) 
ISEED= IX 
ISEED1 = IX 

IF ( . NOT . GOODINPUT) STOP 

OBJVALM = GUESS 

DO 50 1=1 , NUMASSETS 

TLAT { I ) = LATASSET ( I ) 

TLON ( I ) = LONASSET ( I ) 

50 CONTINUE 



RUN THE SIMULATION WITH INITIAL DATA POINTS 

CALL SMUGS IM (TLAT, TLON, TOTCOST) 

OBJVAL0 = TOTCOST 
OBJVALB = OBJVAL0 

IF ( OBJVALB .GT. GUESS) GUESS = 1.5 * OBJVALB 



BEGIN THE ANNEALING 

MOVES = 0 
NOTMOVE = 0 
INFEAS = 0 
FINISH = .FALSE. 

K = 0 

LITTLEN = 0 



100 CONTINUE 



FIND NEW LOCATION FOR ALL ASSETS EXCEPT 
BALLOONS (TYPE 4) 

DO 140 I = 1 , NUMASSETS 

IF ( ASSETTYPE ( I ) . NE . 4) THEN 

120 CONTINUE 

CALL LNORPC ( ISEED1 , A, 2 ) 

DENOM = SQRT { A ( 1 ) * *2 + A(2)**2) 

TLAT(I) = LAT ASSET ( I ) + DELR*A ( 1 ) /DENOM 
TLON { I ) = LONASSET (I) + DELR* A { 2 ) / DENOM 

CHECK FOR FEASIBILITY 

IF ( ASSETLOC ( I ) . EQ . 1) THEN 

CALL LANTCHECK ( TLAT ( I ) , TLON { I ) , FEASIBLE) 
ELSEIF { ASSETLOC { I ) . EQ . 2) THEN 

CALL PACCHECK {TLAT { I ) , TLON { I ) , FEASIBLE) 

ELSE 

CALL SWBCHECK { TLAT ( I ) , TLON ( I ) , FEAS I BLE ) 

ENDIF 

IF ( . NOT . FEASIBLE) INFEAS = INFEAS + 1 
IF ( . NOT . FEASIBLE) GO TO 120 

ENDIF 

140 CONTINUE 

RUN THE SIMULATION WITH NEW DATA POINTS 
CALL SMUGS IM ( TLAT , TLON # TOTCOST ) 

OBJVAL1 = TOTCOST 

DELOBJ = OBJVAL1 - OBJVALO 

ADJUST MAX VALUE GUESS IF REQUIRED 

IF ( OBJVAL1 .GT. OBJVALB) OBJVALB = OBJVAL1 

IF ( OBJVALB .GT. GUESS) GUESS = 1.5 * OBJVALB 

IF (OBJVALB .GT. OBJVALM) OBJVALM = GUESS 

DETERMINE WHETHER OR NOT TO ACCEPT THE MOVE 
IF (DELOBJ .GE. 0.0) THEN 
ACCEPT = .TRUE. 

ELSE 

CK = OBJVALM - OBJVALO 

IF ( CK .EQ. 0.0) CK = 0.000001 

P = EXP (BETA * DELOBJ / CK) 

CALL LRNDPC ( ISEED1 , V, 1 ) 

IF ( P .GE. V) THEN 

ACCEPT = .TRUE. 

ELSE 

ACCEPT = .FALSE. 

K = K + 1 

NOTMOVE = NOTMOVE + 1 

ENDIF 

ENDIF 

IF (ACCEPT) THEN 
COUNT = 0 

DO 160 I = 1 , NUMASSETS 

LAT ASSET ( I ) = TLAT ( I ) 

LONASSET ( I ) = TLON ( I ) 

160 CONTINUE 

MOVES = MOVES + 1 
OBJVALO = OBJVAL1 
LITTLEN = LITTLEN + 1 
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ENDIF 

IF (K .EQ. S) THEN 
DELROLD = DELR 
DELR = DELRLOW 
ELSE IF (K .GE. (2 * S)) THEN 
FINISH = .TRUE. 

ELSEIF ( K .GT. S .AND. ACCEPT) THEN 
DELR = DELROLD 
ENDIF 

IF ( LITTLEN .GE. BIGN) THEN 

OBJVALM = OBJVALM - ALPHA * (OBJVALM - OBJVALB ) 
DELR = DELR - ALPHA * (DELR - DELRMIN) 

ENDIF 

IF (OBJVALM .LT. OBJVALB) OBJVALM = OBJVALB 



I F ( . NOT . FINISH) GO TO 100 
CLOSE (55) 

WRITE THE RESULTS 

OPEN (75, FILE= ' ANNEAL . OUT ' ) 

WRITE(75, 520) OBJVALO/IOOO . 0 
520 FORMAT (IX, 'MAX TOTAL COST TO SMUGGLERS = \F10.2) 
WRITE(75, 530) 

530 FORMAT ( 10X, ' (IN THOUSANDS)') 

WRITE (75, 500 ) MOVES 

500 FORMAT ( IX, 'OPT VALUE FOUND IN ',110,' STEPS') 

WRITE (75, 501) NOTMOVE 

501 FORMAT (I9X, 110, ' STEPS NOT ACCEPTED') 

WRITE (75, *) 

WRITE (75, 510) INFEAS 

510 FORMAT (IX, 110, ' INFEASIBLE MOVES FOUND') 

WRITE (75, * ) 

WRITE (75, * ) 

WRITE (75, * ) 

WRITE (75, 534) 

534 FORMAT ( 1 IX, 'FINAL LOCATIONS') 

WRITE(75, 535) 

535 FORMAT (1 IX, 'ASSET' , ' LATITUDE ', 'LONGITUDE ') 

DO 550 1=1 , NUMASSETS 

WRITE (7 5, 54 0 ) ASSETNAME ( I) , LAT ASSET ( I ) , LONASSET ( I ) 
540 FORMAT (1X,A15,F10.3,F10.3) 

550 CONTINUE 

WRITE (75, *) 

WRITE(75, 589) 

589 FORMAT ( 'FINAL PROBABILITY OF INTERDICTION BY ROUTE') 
WRITE(75, 590) 

590 FORMAT ( 'ROUTE' , 5X, 'ROUTENAME' , 4X, ' PROBINTERDICT ' ) 

DO 600 ROUTENR = 1, NUMROUTES 

WRITE (7 5, 610) ROUTENR, ROUTENAME (ROUTENR) , 

& PROBINTERDICT ( ROUTENR, 1) 

610 FORMAT ( IX , 14 , 5X , A1 0 , 5X, F8 . 4 ) 

600 CONTINUE 
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CPUTIME = DTIME (TARRAY) 

WRITE (75, *) 

WRITE (75, * ) ' CPU TIME = ', CPUTIME, ' SECONDS' 

WRITE (75, * ) ' MAX VALUE FOUND = ' , MAXNOW 
WRITE (75, * ) ' FINAL GUESS VALUE = ' , OBJVALM 
WRITE ( 75 , *) 

WRITE (7 5,*) ' NEXT LOC SEED = ' , ISEED1 

WRITE (75, * ) ' RUN COMPLETED' 

CLOSE (75) 

STOP 

END 

******* SUBROUTINES TO CHECK FEASIBLE REGIONS ******************* 
SUBROUTINE LANTCHECK ( LAT , LON, FEASIBLE) 

★ 

REAL LAT , LON 
LOGICAL FEASIBLE 

★ 

FEASIBLE = .FALSE. 

IF (LON . LE. 87.0 .AND. LON .GE.67.0 .AND. 

& LAT .LE. 25.0 .AND. LAT . GE . 14.0) FEASIBLE = .TRUE. 

RETURN 

END 

★ ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★llr***************-*^*************** 



SUBROUTINE PACCHECK { LAT , LON, FEASIBLE) 

REAL LAT, LON 
LOGICAL FEASIBLE 

FEASIBLE = .FALSE. 



IF (LAT .LE. 0.0) THEN 



ELSE IF ( LON 
IF (LAT , 
ELSEIF (LON 
IF (LAT . 
ELSEIF (LON 
IF (LAT , 
ELSEIF (LON 
IF (LAT . 
ELSEIF (LON 
IF (LAT . 
ELSEIF (LON 
IF (LAT . 
ELSEIF (LON 
IF (LAT . 



& 

END IF 

RETURN 

END 





FALSE . 




.LE. 


120.0 .AND. LON .GE. 


118.691) THEN 


LE. 


33.773) FEASIBLE = .TRUE. 


.GE. 


118.0) THEN 




LE. 


32.50) FEASIBLE = .TRUE. 


.GE. 


115.50) THEN 




LE. 


(-132.70 - 1.40* LON ) ) 


FEASIBLE = .TRUE. 


.GE. 


113.00) THEN 




LE. 


(-64.40 + . 80 *LON ) ) FEASIBLE = .TRUE. 


.GE. 


107.00) THEN 




LE. 


(-69.1667 + .8333 *LON) 


) FEASIBLE = .TRUE 


.GE. 


84.50) THEN 




LE. 


(-39.4444 + . 5 5 5 6 * LON ) 


) FEASIBLE = .TRUE 


.GE. 


79.00) THEN 




LE. 


(-15.5455 + .2727* LON ) 


.AND. 


LAT 


• GE. (121.5 - 1.5* LON ) 


) FEASIBLE = .TRUE 


★ ★ ★ ★ 




■*★★*■■*■*★■*★■*★•*•*■*■*■*■*★ 



r ★ ★ ★ -j 



SUBROUTINE SWBCHECK (LAT, LON, FEASIBLE) 

REAL LAT, LON 
LOGICAL FEASIBLE 
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o o 



FEASIBLE 



.FALSE. 



IF (LON . LE 


. 117 


.13 .AND. LON .GE. 115.0) THEN 




IF (LAT 


.GE. 


(33.74 - .0094* LON ) ) FEASIBLE = . 


TRUE. 


ELSEIF ( LON 


.GE. 


111.0) THEN 




IF (LAT , 


.GE. 


(-7.015 + .345* LON ) ) FEASIBLE = . 


TRUE. 


ELSEIF (LON 


.GE. 


108.295) THEN 




IF (LAT 


.GE. 


(33.7421 - .0222* LON ) ) FEASIBLE = 


.TRUE. 


ELSEIF (LON 


.GE. 


106.37) THEN 




IF (LAT , 


. GE . 


31.742) FEASIBLE = .TRUE. 




ELSEIF { LON 


.GE. 


105.0) THEN 




I F ( LAT . 


.GE. 


(-51.18 + .7796* LON ) ) FEASIBLE = 


.TRUE. 


ELSEIF (LON 


.GE. 


104.474) THEN 




IF (LAT , 


.GE. 


(-173.537 + 1 .9449* LON ) ) FEASIBLE 


= .TRUE 


ELSEIF (LON 


.GE. 


103.474) THEN 




IF (LAT , 


.GE. 


(-38.3616 + .6510* LON ) ) FEASIBLE 


= .TRUE. 


ELSEIF (LON 


.GE. 


103.36) THEN 




IF (LAT . 


.GE. 


29.0) FEASIBLE = .TRUE. 




ELSEIF (LON 


.GE. 


102. 731) THEN 




IF (LAT , 


.GE. 


(166. 5411 - I . 3 3 07* LON) ) FEASIBLE 


= .TRUE 


ELSEIF ( LON 


.GE. 


101.462) THEN 




IF (LAT . 


.GE. 


29.837) FEASIBLE = .TRUE. 




ELSEIF (LON 


.GE. 


100.747) THEN 




IF (LAT . 


.GE. 


(-72.6813 + I . 0098*LON) ) FEASIBLE 


= .TRUE 



ELSE IF ( LON 
I F ( L AT , 
ELSE IF (LON 
I F { L AT . 
ELSE IF ( LON 
I F ( LAT . 
ENDIF 
RETURN 
END 



.GE. 99.64) THEN 

GE. (-104.1222 + 1 .3225* LON ) ) FEASIBLE = .TRUE. 
.GE. 99.00) THEN 

GE. (-168.6709 + 1.9703*LON)) FEASIBLE = .TRUE. 
.GE. 97.133) THEN 

GE. (1.7858 + .2485* LON ) ) FEASIBLE = .TRUE. 






SMUGS IM SUBROUT I NE 






SUBROUT I NE SMUGS IM ( TLAT , TLON , TOTCOST ) 

C SMUGGLERS ' SIMULATION - 2/3/87 

* ADAPTED FROM RAND CORPORATION SIMULA! L U <4 aD/UU i 1 

* RESPONSE : A MODEL OF DRUG INTERDICTION . [Ref. 2] 

C INSTALLATION REMARKS: 

C 

C YOU'LL ALSO NEED TO SELECT "SIZING PARAMETERS" LARGE ENOUGH 

C TO HANDLE YOUR ANALYSIS , BUT NOT SO LARGE AS TO MAKE THE 

C PROGRAM TOO LARGE TO RUN ON YOUR MACHINE. SEE THE BEGINNING 

C OF COMMON. FOR FOR PARAMETERS. 

£★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 
C PROGRAM VARIABLES: 



INCLUDE 'COMMON. FOR' 

REAL TLAT ( MAXASSETS ) , TLON (MAX ASSETS ) 
INTEGER* 4 T, R 

INITIALIZE THE DATA COLLECTION ARRAYS. 
DO 30, R= 1 , MAXROUTES 
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AMOUNT ATTEMPTED ( R ) = 0. 
AMOUNTSUCCEEDED(R) = 0. 

ATTEMPTS (R) = 0. 

SUCCESSES (R) = 0. 

DO 20, THISPHASE = 1 , MAXPHASES 

ATTEMPTSBY PHASE ( R , TH IS PHASE) = 0. 
SUCCESSESBYPHASE(R, THISPHASE) = 0. 
CONTINUE 



30 CONTINUE 

SUCCESSCOSTS=0 . 0 
FAILURECOSTS=0 . 0 

SET UP THE TABLE OF WEIGHTS FOR PAST SHIPMENTS. 

KMEMORY = ALOG { MEMORY VALUE ) / { -LONG PAST) 

DO 40, T=0, LONG PAST 

EX PT ABLE (T ) = EXP { -KMEMORY * REAL { T ) ) 

40 CONTINUE 

CALL THE SUBROUTINE TO COMPUTE ROUTE PROBABILITIES BASED ON 

ASSET LOCATIONS 

CALL PROBCOMP ( TLAT , TLON ) 

CALL THE SUBROUTINE THAT DOES THE RUN IN. 

CALL PEACE 

DO EACH TRIAL. 

DO 50, NTRIAL=1 , NUMTRIALS 
CALL SIMULATE 
50 CONTINUE 

COMPUTE OUTPUT TO ANNEALING ALGORITHM 

INITIALIZE THE TOTALS. 

TOTATTEMPTED=0 . 

TOTSH I PPED=0 . 

TOTATTEMPTS =0 . 

TOTSUCCESSES=0 . 

TOTSUCCESSCOST=0 . 

TOTFAI LURECOST=0 . 



& 

Sc 

& 



30 



COMPUTE THE TOTAL NUMBER OF ATTEMPTS, SUCCESSES, COSTS, ETC. 
DO 30, R=1 , NUMROUTES 

TOTATTEMPTS = TOTATTEMPTS + ATTEMPTS (R) 

TOTSUCCESSES = TOTSUCCESSES + SUCCESSES (R) 

TOTATTEMPTED = TOTATTEMPTED + 

AMOUNTATTEMPTED { R ) 

TOTSHIPPED = TOTSHIPPED + AMOUNTSUCCEEDED ( R ) 

SUCCESSCOSTS = SUCCESSCOSTS + 

AMOUNTSUCCEEDED (R) * DRUGCOST 

FAI LU RECOSTS = FAILURECOSTS + 

{ AMOUNTATTEMPTED (R) - 

AMOUNTSUCCEEDED (R) ) * 

DRUGCOST 

CONTINUE 

TOTSUCCESSCOST=TOTSUCCESSCOST+ SUCCESSCOSTS 
TOTFA I LURECOST=TOTFAILURECOST+ FAILURECOSTS 



C 

C 



GET THE AVERAGE NUMBER OF ATTEMPTS, SUCCESSES, COSTS, ETC. 
PER TRIAL. 



53 



non ooooo o nonoo 



TOTATTEMPTS = TOTATTEMPTS/NUMTRI ALS 
TOTSUCCESSES = TOTSUCCESSES/NUMTRIALS 
TOTSUCCESSCOST = TOTSUCCESSCOST/NUMTRI ALS 
TOTFAILURECOST = TOTFAILURECOST/NUMTRI ALS 
TOTCOSTS = TOTSUCCESSCOST + TOTFAILURECOST 



RETURN 

END 

£★**★★*★*****★★*****★**★********★*★*★*****★***★*******★★★★★*★*★★****★ 

SUBROUTINE COMPFACT (RNUMERATOR, RDENOMINATOR, R) 

GET THE NUMERATOR AND DENOMINATOR OF THE ”R" FACTOR, WHICH 
WILL BE USED TO INCREASE THE PROBABILITY OF INTERDICTION 
ON ROUTES WITH HIGHER THAN AVERAGE TRAVEL. 

★ ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★•A-**** 

PROGRAM VARIABLES: 

INCLUDE ' COMMON . FOR ' 

INTEGER* 4 R , RNUMERATOR , RDENOMINATOR, S 

★ ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★-a-* 

RNUMERATOR = 0 
RDENOMINATOR = 0 
DO 10, S=0 , LONGPAST 

RDENOMINATOR = RDENOMINATOR + PASTSHIPMENTS (R, S) 

10 CONTINUE 

DO 20, S=0, RECENT PAST 

RNUMERATOR = RNUMERATOR+ PASTSH1 PMENTS ( R , S ) 

20 CONTINUE 
RETURN 
END 



★ ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★•A** 



SUBROUTINE GETDATA (GOODINPUT) 

READ IN THE INPUT DATA. 

PROGRAM VARIABLES: 

INCLUDE ' COMMON . FOR ' 

REAL STEPSIZEM 

I NTEGER * 4 DAY 1 , DAY 2 , DAYN , M , N , TEMP , TEMP2 , 
& P, ROUTENR , I, J, PTS {MAXPTS ) , 

Zc METHOD, 

Zc POINT, TYPE, ASSET 

CHARACTER* 15 PO I NTNAME {MAX POINTS ) 
LOGICAL GOODINPUT 

LOGICAL TEMPIF 



GOODINPUT WILL INDICATE WHETHER THE DATA WAS CLEAN AND THE 
SIMULATION SHOULD BE RUN. 

GOODINPUT = .TRUE. 
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OPEN ( 55 , FILE= ' ANNEAL . DAT ' ) 

★ 

READ ( 55 , *) 

READ (55, *) 

READ ( 55 , 4 010 ) STEPS I Z EM 
READ ( 55 , 4010} STEPS I ZEMIN 
READ(55, 4010) STEPSIZELOW 
4 010 FORMAT ( 2 OX , F5 . 1 ) 

* APPROXIMATION OF STEP SIZE IN DEGREES 
DELR = STEPSIZEM/60 . 0 

DELRMIN = STEPSIZEMIN/60 . 0 
DELRLOW = STEPSIZELOW/ 60.0 

★ 

READ(55, 4020) BETA 
4 02 0 FORMAT (2 OX, F5. 2) 

★ 

READ ( 55 , 4 03 0 ) S 
4030 FORMAT (2 OX, 15) 

★ 

READ(55, 4040) ALPHA 
4040 FORMAT ( 20X, F5. 3) 

★ 

READ(55, 4050) GUESS 
4050 FORMAT ( 2 OX , FI 5 . 2 ) 

★ 

READ(55, 4060) BIGN 
4060 FORMAT (2 OX, 110) 

★ 

CLOSE (55) 



C READ IN THE OVERALL SOAR SIMULATION DATA. 

OPEN (7 6, FILE= 'INPUT. DAT ' ) 

READ (76, *) 

READ (76,5010) NUMMETHODS , NUM PHASES , NUMROUTES , 

& NUMTRIALS, RUN IN, ENDTIME, MEMORYVALUE, IX 

5010 FORMAT (6 (15, 4X) , F10.5, 2X, 15) 

C CHECK WHETHER TOO MANY METHODS HAVE BEEN REQUESTED. 

TEMPI F = (NUMMETHODS .GT. MAXMETHODS) 

CALL WRTERROR (TEMPIF, GOODINPUT) 

C CHECK WHETHER TOO MANY PHASES HAVE BEEN REQUESTED. 

TEMPIF = (NUMPHASES .GT. MAXPHASES) 

CALL WRTERROR (TEMP IF, GOODINPUT) 

C CHECK WHETHER TOO MANY ROUTES HAVE BEEN REQUESTED. 

TEMPIF = (NUMROUTES .GT. MAXROUTES) 

CALL WRTERROR (TEMPIF, GOODINPUT) 

C CHECK WHETHER TOO MANY DAYS HAVE BEEN REQUESTED. 

TEMPIF = (ENDTIME .GT. MAXDAYS ) 

CALL WRTERROR (TEMPIF, GOODINPUT) 

C MAKE SURE THE INITIAL SEED IS ODD. (IN CASE WE RUN ON A SUN.) 

IX = IX / 2 
IX = IX * 2 + 1 

C READ IN THE DRUG RELATED DATA RECORDS. 

READ (76, * ) 
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READ (76,5020) DRUGNAME, DRUGCOST, EXSHI PMENTINTRVL, 

£ EXSHI PMENTSIZE, DAILYAMOUNT 

5020 FORMAT (A10.4F10.5) 

C READ IN THE METHOD RELATED RECORDS. 

TEMP = NUMMETHODS 

IF (TEMP .GT. MAXMETHODS) TEMP = MAXMETHODS 
READ (7 6, *) 

DO 40, M= 1 , TEMP 

READ (76,5030) METHODNAME (M) , 

£ RISKCOMP (M ) , RISKCOMPEXP(M) ,CAPCOST(M) , 

£ CAPACITY (M) 

5030 FORMAT (A5, 5X, 4F10 . 5) 

40 CONTINUE 

IF (NUMMETHODS .GT. TEMP) THEN 
DO 50, M=TEMP+1, NUMMETHODS 
READ (76,5030) 

50 CONTINUE 

ENDIF 

C READ IN THE DAYS WHEN EACH PHASE ENDS. SET UP THE VECTOR 

C INDICATING WHICH PHASE IS IN EFFECT FOR EACH DAY. 

DAY1 = 0 

TEMP = NUMPHASES 

IF (TEMP .GT. MAXPHASES) TEMP = MAXPHASES 
DO 70, N=1 , TEMP 

READ (76,5040) DAY 2 
5040 FORMAT (15) 

TEMP2 = DAY2 

IF (TEMP2 .GT. MAXDAYS) TEMP2 = MAXDAYS 
IF (TEMP2 .GE. DAY1 ) THEN 
WRITE (77,6016) N, TEMP2 

6016 FORMAT (8X,‘ PHASE ',12,' LASTS THROUGH DAY ',14) 

DO 60, DAYN=DAY1 , TEMP2 

CURRENTPHASE ( DAYN) = N 
60 CONTINUE 

DAY1 = TEMP2 + 1 
ENDIF 

70 CONTINUE 
CLOSE (76) 

* READ IN THE NETWORK ROUTES AND ROUTE DATA 

OPEN ( 96 , FILE= ' ROUTES . DAT ' ) 

READ (96, *) 

READ (96, *) 

DO 250 ROUTENR = l.NUMROUTES 
DO 210 P = 1, MAXPOINTS 
NEXT (ROUTENR, P) = 0 
210 CONTINUE 

READ (96,5075) ROUTENAME ( ROUTENR ) , ROUTECOST ( ROUTENR ) , 

£ ROUTEMETHOD ( ROUTENR ) , ROUTELOC ( ROUTENR ) 

507 5 FORMAT (9X, A10, 3X, F10.2, 2X, 15, 15) 

DO 220 1=1 , MAXPTS 
PTS(I) =0 
220 CONTINUE 

READ ( 96 , 5080 ) FIRST ( ROUTENR ) , ( PTS ( I ) , 1=1, MAXPTS) 

5080 FORMAT (10X, 915) 

P=FIRST (ROUTENR) 
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J = 1 

IF ( P .NE. 0) THEN 

NEXT (ROUTENR, P) = PTStJj 
P = NEXT ( ROUTENR, P) 

J = J + 1 
GO TO 240 
ENDIF 
250 CONTINUE 
CLOSE (96) 

* READ IN THE ASSET DATA 

OPEN ( 98, FILE=' ASSET. DAT' ) 

READ ( 98 , 5300) NUMASSETS 
READ (98, 5300 ) NUMASSETTYPES 
5300 FORMAT ( 2 5X , 15) 

C CHECK WHETHER TOO MANY ASSETS HAVE BEEN REQUESTED. 

TEMP IF = (NUMASSETS .GT. MAXASSETS) 

CALL WRTERROR (TEMPIF, GOODINPUT) 

C CHECK WHETHER TOO MANY ASSET TYPES HAVE BEEN REQUESTED. 

TEMPIF = (NUMASSETTYPES .GT. MAXASSETTYPES ) 

CALL WRTERROR (TEMPIF, GOODINPUT) 

READ (98, *) 

DO 350 ASSET = 1, NUMASSETS 

READ ( 98 , 532 0) ASSETNAME (ASSET ) , 

£ LAT ASSET ( ASSET ) , LONASSET ( ASSET ) , 

& ASSETTY PE ( ASSET ) , ASSETLOC ( ASSET ) 

5320 FORMAT (5X, A15,2(F10.3),2(I5)) 

350 CONTINUE 
READ (98, *) 

READ (98, *) 

READ ( 98 , *) 

READ (98, *) 

DO 360 TYPE=1 , NUMASSETTYPES 

READ(98, 5325) ( RANGE (TYPE, METHOD) , METHOD=l , NUMMETHODS) 

5325 FORMAT (15X, 3 (F5. 1, 5X) ) 

360 CONTINUE 
READ (98, * ) 

READ (98, *) 

READ (98, *) 

READ (98, *) 

DO 390 TYPE=1 , 4 

READ(98, 5330) ( IPD(TYPE, METHOD) ,METHOD=l, 3) 

5330 FORMAT ( 1 5X , 3 ( F5 . 3 , 5X ) ) 

390 CONTINUE 
CLOSE (98) 

* READ IN THE GEOGRAPHICAL DATA 

OPEN ( 95, FI LE=' GEOGR.DAT' ) 

READ (95, *) 

392 READ ( 95 , 5400 , END=395 ) POINT, POINTNAME ( POINT ) , LAT ( POINT ) , LON(POINT) 
5400 FORMAT ( IX, 1 5, 5X, A15 , 5X, F7 . 3 , 5X, F7 . 3 ) 

GO TO 392 
395 CONTINUE 
CLOSE (95) 

RETURN 

END 
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SUBROUTINE INITSIM 

C INITIALIZE FOR THE CURRENT SAMPLE POINT. 



C PROGRAM VARIABLES: 

INCLUDE ' COMMON . FOR ' 
INTEGER* 4 D, R 
REAL RANDNM 



C SET THE SHIPMENT COUNTER TO ZERO. 

NUMSHIPMENT = 0 



C DETERMINE WHEN THE NEXT SHIPMENT WILL OCCUR. 

CALL LRNDPC ( IX, RANDNM, 1 ) 

NEXTEVENT = ENDTIME +1.0 

NEXTSHIPMENT = - EXSH I PMENT INTRVL * ALOG ( RANDNM) 
IF (NEXTSHIPMENT . LT . NEXTEVENT) THEN 
NEXTEVENT = NEXTSHIPMENT 
ENDIF 



C GET THE PAST FROM THE LONG RUN IN. 

DO 30, D= 0 , LONGPAST 

DO 20, R=1 , NUMROUTES 

PASTSHIPMENTS (R, D) = PEACESHI PMENTS ( R , D) 
PASTFAILURES ( R , D) = PEACEFAILURES (R, D) 

20 CONTINUE 

30 CONTINUE 
RETURN 
END 









SUBROUTINE PEACE 
C DO THE INITIAL RUN IN. 

£★**★**★****★★★***★*★*★★*★*★*★★*★**★*★****★***★★★**★★★•*★★*★*★★**★★*■*★■*** 
C PROGRAM VARIABLES: 



INCLUDE ' COMMON . FOR ' 

REAL AMOUNT, PROBCAUGHT, RFACTOR, TEMPTIME, RANDNM, ACTUALRISKCOMP 
INTEGER *4 D, R, RDENOMINATOR , ENUMERATOR, T , TRI PS 
LOGICAL TRIPSUCCESS 






C DETERMINE THE TIME OF THE FIRST SHIPMENT. 

NEXTEVENT = RUNIN +1.0 
CALL LRNDPC (IX, RANDNM, 1) 

NEXTSHIPMENT = -EXSH I PMENT INTRVL * ALOG (RANDNM) 
IF (NEXTSHIPMENT . LT . NEXTEVENT) THEN 
NEXTEVENT = NEXTSHIPMENT 
ENDIF 
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220 

230 



INITIALIZE THE ARRAYS DESCRIBING THE PAST. 
DO 230 , D=0 , LONGPAST 

DO 220 , R=1 , NUMROUTES 

PASTSH I PMENTS { R , D ) = 0 
PASTFAILURES (R, D) = 0 
CONTINUE 
CONTINUE 



FOR EACH DAY OF RUNIN.. 
DO 100, DAYNOW=l, RUNIN 
THISPHASE = 1 



10 



20 



SHIFT THE ARRAYS DESCRIBING THE PAST. 

DO 20, R=l, NUMROUTES 

DO 10, D=LONGPAST, 1,-1 

PASTSH I PMENTS (R, D) = PASTSH I PMENTS { R, D-l ) 
PASTFAILURES (R,D) = PASTFAILURES (R, D-l ) 
CONTINUE 

PASTSH I PMENTS {R, 0) = 0 
PASTFAILURES ( R, 0 ) = 0 
CONTINUE 



30 



LOOP THROUGH THE DAYS SHIPMENTS. 

IF (NEXTEVENT .GE. DAYNOW + 1.0) GO TO 100 



GET THE AMOUNT OF THE NEXT SHIPMENT. 

AMOUNT = EXSH I PMENTS I ZE 

SELECT THE ROUTE TO BE USED. 

CALL SELROUTE (R, TRIPS, AMOUNT, ACTUALR I SKCOMP) 

COMPUTE THE NUMERATOR AND DENOMINATOR OF THE 
" R " FACTOR. 

CALL COMPFACT ( RNUMERATOR , RDENOMINATOR , R ) 

FOR EACH TRIP REQUIRED TO GET AMOUNT SHIPPED. . . 

DO 70, T=l, TRIPS 

COMPUTE THE "R" FACTOR. 

RFACTOR =1.0 

IF {RDENOMINATOR .GT. 0.0) THEN 

RFACTOR = RNUMERATOR /RDENOMINATOR 
TEMPTIME = NEXTEVENT 
IF {TEMPTIME .GT. RECENTPAST) THEN 
IF (TEMPTIME .GT. LONGPAST) 

TEMPTIME = LONGPAST 

RFACTOR = RFACTOR * TEMPT IME/RECENTPAST 
END IF 
END IF 

IF (RFACTOR . LT . 1.0) RFACTOR =1.0 

COMPUTE THE PROBABILITY OF INTERDICTION. 

IF (PROBINTERDICT (R, THISPHASE) .GE. .9999) THEN 
PROBCAUGHT =1.0 

ELSE 

PROBCAUGHT = 1.0 - EXP ( RFACTOR* ALOG ( 1 . 0 - 
PROBINTERDICT (R, THISPHASE) ) ) 

ENDIF 



DETERMINE WHETHER THE TRIP WAS SUCCESSFUL . 
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CALL LRNDPC ( IX , RANDNM, 1 ) 

TRIPSUCCESS = (RANDNM .GE. PROBCAUGHT) 



1 

+ 1 



C GET THE TIME AND TYPE OF THE NEXT SHIPMENT. 

CALL LRNDPC ( IX , RANDNM, 1 ) 

NEXTSHIPMENT = NEXTEVENT - 

Sc EXSHI PMENTINTRVL * ALOG ( RANDNM) 

NEXTEVENT = RUNIN + 1.0 
IF (NEXTSHIPMENT . LT . NEXTEVENT) THEN 
NEXTEVENT = NEXTSHIPMENT 
ENDIF 
GO TO 30 
100 CONTINUE 



C DO THE BOOKKEEPING. (NOT VERY EXTENSIVE 

C DURING THE RUN IN.) 

PASTSH I PMENTS ( R , 0 ) = PASTSHI PMENTS ( R , 0 ) + 
IF (.NOT. TRIPSUCCESS) 

Sc P ASTFAI LURES (R, 0) = PASTFAILURES ( R , 0 ) 

R NUMERATOR = RNUMERATOR + 1 
RDENOMINATOR = RDENOMINATOR + 1 
70 CONTINUE 



C SAVE THE LAST LONGPAST DAYS FOR USE INITIALIZING EACH TRIAL. 

DO 130, D=0, LONGPAST 

DO 120, R=1 , NUMROUTES 

PEACESH I PMENTS (R, D) = PASTSHI PMENTS ( R, D) 

PEACEFAILURES (R, D) = PASTFAI LURES ( R , D) 

120 CONTINUE 

130 CONTINUE 
RETURN 
END 

£★*★*★★*★**★**★★★*★****★**★***★*★★****★★★★*★★★★★★★★★**★★★★★★*★*★★★**★ 






SUBROUT I IME SELROUTE ( RCHOSEN , TR I PS , AMOUNT , ACTUALR ISKCOMP) 

C SELECT THE ROUTE FOR THE NEXT SHIPMENT. 

C PROGRAM VARIABLES: 

INCLUDE ' COMMON . FOR ' 

REAL AMOUNT , CUMPROB , ROUTEPROB, COSTROUTE (MAXROUTES ) , 

Sc ACTUALR ISKCOMP , PROBCAUGHT ( MAXROUTES ) , 

Sc TEMPNUMERATOR (MAXROUTES) , TEMPDENOMINATOR (MAXROUTES ) , 

Sc TOTCOST , WEIGHT , TEMPRISKCOMP (MAXROUTES ) 

INTEGER* 4 RCHOSEN, TRI PS , RMETHOD, R, S 






C INITIALIZE THE NUMERATOR AND DENOMINATOR FOR EACH ROUTE. 

DO 10, R=l, NUMROUTES 

TEMPNUMERATOR (R) =0.0 
TEMPDENOMINATOR (R) =0.0 
COSTROUTE (R) =0.0 
10 CONTINUE 

C FOR EACH DAY TO BE CONSIDERED, ADD IN ITS CONTRIBUTION. 
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DO 30, S=0, LONG PAST 

WEIGHT = EXPTABLE ( S ) 

DO 20, R=1 , NUMROUTES 

TEMPDENOMINATOR (R) = 

& TEMPDENOMINATOR (R) + WEIGHT* PASTSHI PMENTS (R, S) 

TEMPNUMERATOR ( R ) = TEMPNUMERATOR ( R ) + 

& WEIGHT* PASTFAILURES (R, S) 

20 CONTINUE 
30 CONTINUE 



FOR EACH ROUTE, COMPUTE THE PERCEIVED PROBABILITY OF BEING 
CAPTURED AND HENCE THE EXPECTED COST OF USING THE ROUTE. 

THE PROBABILITY A ROUTE WILL BE CHOSE WILL BE PROPORTIONAL 
TO THE INVERSE EXPECTED COST OF USING THE ROUTE, SO GET THE 
TOTAL OF THE INVERSE EXPECTED COSTS OF USING EACH ROUTE. 
TOTCOST = 0.0 
DO 40, R=1 , NUMROUTES 

RMETHOD = ROUTEMETHOD ( R ) 

TRIPS = INT (AMOUNT/CAPACITY (RMETHOD) + 0.999) 

IF (TEMPDENOMINATOR (R) .GT. 0.001) THEN 

PROBCAUGHT (R) = TEMPNUMERATOR ( R ) /TEMPDENOMINATOR ( R ) 

ELSE 

PROBCAUGHT (R) =0.0 
END IF 

IF (RISKCOMPEXP( RMETHOD) . LE . 0.0001) THEN 
TEMPR ISKCOMP ( R ) =R I SKCOMP ( RMETHOD) 

ELSE 

TEMPRISKCOMP (R) =RISKCOMP ( RMETHOD) * 

£c ( (2 ‘PROBCAUGHT (R) ) * ‘RISKCOMPEXP ( RMETHOD) ) 

ENDIF 

COSTROUTE (R) = 1 . 0/ (TRIPS* ( PROBCAUGHT (R) *CAPCOST (RMETHOD) + 
£< ROUTECOST (R) + TEMPRISKCOMP ( R )) + 

U PROBCAUGHT ( R ) * AMOUNT * DRUGCOST ) 

TOTCOST = TOTCOST + COSTROUTE (R) 

40 CONTINUE 

NOW CHOOSE THE ROUTE, WHERE THE PROBABILITY OF CHOOSING A 
GIVEN ROUTE IS PROPORTIONAL TO THE INVERSE COST OF USING 
THAT ROUTE. 

CALL LRNDPC ( IX, ROUTEPROB, 1 ) 

CUMPROB = 0.0 
R = 0 

50 CONTINUE 

R = R + 1 

CUMPROB = CUMPROB + COSTROUTE ( R ) /TOTCOST 
IF (CUMPROB . LT. ROUTEPROB) GO TO 50 
RCHOSEN = R 

RMETHOD = ROUTEMETHOD (RCHOSEN) 

TRIPS = INT (AMOUNT/CAPACITY (RMETHOD) + 0.999) 
ACTUALRISKCOMP=TEMPRISKCOMP ( RCHOSEN) 



RETURN 

END 






*★*★★★★*★★★★★*★★★★★★★*★★*★*★★*★★★*★**★*★***★**★★**********★★*★*★**★★ 



SUBROUTINE SIMULATE 

THIS SUBROUTINE PERFORMS THE SIMULATION FOR A SINGLE SAMPLE 
POINT. 
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C PROGRAM VARIABLES: 



INCLUDE 'COMMON. FOR' 

REAL AMOUNT, PROBCAUGHT, RFACTOR, TEMPTIME, RANDNM, ACTUALRI SKCOMP , 
& T IMESH I PPED ( MAXSH I PMENTS ) 

INTEGER* 4 D, R , RDENOMINATOR, RNUMERATOR, T, TRIPS, 
it ROUTEUSED (MAXSH I PMENTS ) 

LOGICAL TRI PSUCCESS , SUCCESS (MAXSHI PMENTS ) 



£★**★★*★★★★■*★★★★★★★*★*★*★★★★★■**★★★★★★•*★*★★★*★★★■*■*■*•*★■*★•*★★■*•****■**•*■*•**■*■■*■** 

C INITIALIZE FOR THIS TRIAL. 

CALL INITSIM 



C FOR EACH DAY . . . 

DO 100, DAY NOW =0, ENDTIME 

C GET THE POINTER INTO THE ARRAY OF INTERDICTION PROBABILITIES. 

THISPHASE = CURRENT PH ASE(DAYNOW) 

C SHIFT THE PAST HISTORY ARRAYS. 

DO 20, R=1 , NUMROUTES 

DO 10, D= LONG P AST, 1, -1 

PASTSHI PMENTS ( R, D) = PASTSH I PMENTS ( R , D- 1 ) 
PASTFAILURES {R, D) = PASTFAILURES (R, D-l ) 



10 


CONTINUE 








PASTSH I PMENTS (R, 0) : 


= 0 






PASTFAILURES (R, 0) = 


0 




20 


CONTINUE 








LOOP THROUGH THE EVENTS 


THAT HAPPEN 


TODAY . 


30 


IF (NEXTEVENT .GE. DAYNOW + 1.0) GO 


TO 10 0 



DETERMINE THE NEXT TYPE OF DRUG TO BE SHIPPED AND 
HOW MUCH. 

AMOUNT = EXSH I PMENTS I ZE 

C SELECT THE ROUTE. 

CALL SELROUTE ( R , TR I PS , AMOUNT , ACTU ALR I SKCOM P ) 

C GET THE NUMERATOR AND DENOMINATOR OF THE "R" FACTOR. 

CALL COMPFACT ( RNUMERATOR, RDENOMINATOR , R ) 

C MULTIPLE TRIPS MAY BE REQUIRED TO SHIP THE GIVEN 

C AMOUNT OF DRUG ON THE SELECTED ROUTE, DUE TO 

C CAPACITY CONSTRAINTS. THE VARIABLE TRIPS CONTAINS 

C THE NUMBER OF TRIPS THAT WILL BE REQUIRED. I EXPECT 

C THAT THIS VARIABLE WILL USUALLY BE CONTAIN A ONE 

C AND THAT AS A RESULT DO-LOOP 70 WILL USUALLY BE 

C EXECUTED ONLY ONCE. 

DO 70, T=l, TRIPS 

C COMPUTE THE "R" FACTOR. 

RFACTOR =1.0 

IF (RDENOMINATOR .GT. 0.0) THEN 

RFACTOR = RNUMERATOR /RDENOMINATOR 
TEMPTIME = NEXTEVENT 
IF (TEMPTIME .GT. RECENTPAST) THEN 
IF (TEMPTIME .GT. LONG PAST) 
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TEMPTIME = LONG PAST 

RFACTOR = RFACTOR * TEMPTIME/ RECENTPAST 
ENDIF 
END IF 

IF {RFACTOR . LT . 1.0) RFACTOR =1.0 

COMPUTE THE PROBABILITY OF INTERDICTION. 

COMPUTE THE PROBABILITY OF INTERDICTION. 

IF ( PROBINTERDICT ( R, THISPHASE) .GE. .9999) THEN 
PROBCAUGHT =1.0 

ELSE 

PROBCAUGHT = 1.0 - EX P ( RFACTOR * ALOG (1.0 - 
PROBINTERDICT { R, THISPHASE ) ) ) 

ENDIF 

DETERMINE WHETHER THE TRIP WAS SUCCESSFUL . 

CALL LRNDPC {IX, RANDNM, 1 ) 

TRIPSUCCESS = (RANDNM .GE. PROBCAUGHT) 

DO THE REQUIRED BOOKKEEPING. 

NUMSHIPMENT = NUMSHIPMENT + 1 
SUCCESS { NUMSHI PMENT ) = TRIPSUCCESS 
T I MESH I PPED { NUMSHIPMENT ) = NEXT EVENT 
ROUTEUSED (NUMSHIPMENT ) = R 
ATTEMPTS (R) = ATTEMPTS (R) +1.0 
AMOUNT ATTEMPTED { R ) = 

AMOUNT ATTEMPTED { R ) + AMOUNT /TRIPS 
ATTEMPTSBYPHASE ( R , THISPHASE) = 

ATTEMPTSBY PHASE { R , THISPHASE ) + 1 
IF (TRIPSUCCESS) THEN 

SUCCESSCOSTS=SUCCESSCOSTS+ ACTUALR ISKCOMP+ 

ROUTECOST (R) 

SUCCESSES ( R ) = SUCCESSES ( R ) + 1.0 
AMOUNTSUCCEEDED (R) = 

AMOUNTSUCCEEDED ( R ) + AMOUNT /TRIPS 
SUCCESSESBYPHASE (R, THISPHASE) = 

SUCCESSESBYPHASE { R , TH ISPHASE) + 1 

ELSE 

PASTFAILURES ( R, 0 ) = PASTFAILURES ( R, 0 ) + 1 

FAILURECOSTS=FAILURECOSTS+ACTUALRISKCOMP+ 
ROUTECOST (R) + 

CAPCOST ( ROUTEMETHOD { R ) ) 

ENDIF 

PASTSH I PMENTS { R , 0 ) = PASTSH I PMENTS { R , 0 ) + 1 
RNUMERATOR = RNUMERATOR + 1 
R DENOMINATOR = R DENOMINATOR + 1 
CONTINUE 

DETERMINE WHEN THE NEXT SHIPMENT OF THIS DRUG WILL BE. 
CALL LRNDPC ( IX , RANDNM , 1 ) 

NEXTSHIPMENT = NEXTEVENT - 

EXSHI PMENTINTRVL * ALOG (RANDNM) 

DETERMINE WHAT THE NEXT SHIPMENT WILL BE. 

NEXTEVENT = ENDTIME + 1.0 
IF (NEXTSHIPMENT . LT . NEXTEVENT) THEN 
NEXTEVENT = NEXTSHIPMENT 
ENDIF 
GO TO 30 
CONTINUE 



u * 



RETURN 

END 

£★★★★*★*★★★★★★★★*★★★★★★*★****★★*★★★★★★★*★**★*★★★**★****■****★*★*★*★★★★ 

£★★*★★********★★★**★★★*★*★*★★*★★★*****★★*★★★★*★★★★★★★★★★*★*★**★*★★*★* 



SUBROUTINE WRTERROR ( ERRORCONDITION , GOODINPUT) 

THIS SUBROUTINE WRITES THE ERROR MESSAGE, IF NEEDED. 

^★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★***********^************** 
C PROGRAM VARIABLES: 

LOGICAL ERRORCONDITION, GOODINPUT 

£★*★*★****★***★★★**★*★★*★*★*★★*★******★★★**★*★*★★***★*★*★★****★**★*★★★★■* 



IF (ERRORCONDITION) THEN 
WRITE (77,6001) 

60 01 FORMAT (' * * * ERROR* * * ' ) 

GOODINPUT = .FALSE. 

END IF 

RETURN 

END 






★★★*★★*★*★★*★***★*★*★★★★★*★★★****★★*★★★★★*★**★★★*★******★★★*★★*★*★★★★ 
SUBROUTINE PROBCOMP ( TLAT , TLON ) 

* THIS SUBROUTINE COMPUTES PROBABILITIES OF INTERDICTION ON 

* ROUTES AS A FUNCTION OF ASSET DISTANCE FROM ROUTE. 

* PROBABILITY OF INTERDICTION IS INPUT TO SMUGSIM. 

* 

INCLUDE ' COMMON . FOR ' 

★ 

REAL Ml , M2 , B1 , B2 , TEMPPI , LATINT, LONINT, DST, DEND, LENGTH, DIST, A, 

& TLAT (MAXASSETS) , TLON (MAX ASSETS ) ,GSDIS,FCTR 

INTEGER P, Q, R 



* COMPUTE PROBABILITY OF INTERDICTION FOR EACH ROUTE 

DO 405 R=1 , NUMROUTES 
DO 400 N=1 , NUMPHASES 

PROB INTERDICT (R, N) =0.0 
400 CONTINUE 

405 CONTINUE 

* FOR EACH ROUTE COMPUTE EQUATIONS FOR ROUTE SEGMENTS 

DO 440 R=1 , NUMROUTES 
P=FIRST ( R ) 

Q=NEXT ( R, P) 

410 CONTINUE 

IF ( LON ( P) .NE. LON (Q) .AND. 

& LAT(P) .NE. LAT (Q) ) THEN 

Ml = ( LAT ( P) - LAT (Q) )/ 

& ( LON ( P ) - LON (Q) ) 

M2=-l .0/Ml 

B1=LAT(P) - Ml *LON ( P) 

END IF 
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TEMPPI=0 . 0 

DO 420 ASSET=1 , NUMASSETS 

IF (LON (P) .NE. LON(Q) .AND. 

Sc LAT(P) .NE. LAT (Q) ) THEN 

B2=LAT ASSET (ASSET) - M2 * LONASSET ( ASSET ) 
LATINT= ( B1 *M2 - B2 *M1 ) / (M2 -Ml ) 

LON I NT = ( B1 -B2 ) / (M2 -Ml ) 

ELSE IF ( LON ( P ) . EQ . LON (Q) ) THEN 

LATINT=LATASSET ( ASSET ) 

LONINT=LON ( P ) 

ELSE 

LATINT=LAT ( P ) 

LONINT=LONASSET (ASSET) 

END IF 



* DETERMINE CLOSEST POINT BETWEEN ASSET AND ROUTE SEGMENT 
DST=GSDIS ( LON { P ) , LAT(P) , 

& LON I NT , LAT I NT ) 

DEND=GSDIS ( LON (Q) , LAT (Q) , 

Sc LON I NT, LAT I NT) 

LENGTH =GSD IS ( LON ( P) , LAT(P) , 

Sc LON(Q) , LAT(Q) ) 

I F ( DST .LE. LENGTH .AND. DEND . LE . LENGTH) THEN 
DIST=GSDIS ( LONASSET (ASSET ) , LATASSET ( ASSET ) , 

Sc LONINT, LATINT) 

ELSEIF ( DST . LT . DEND) THEN 

DIST=GSDIS (LONASSET (ASSET) , LATASSET ( ASSET ) , 

Sc LON { P ) , LAT { P ) ) 

ELSE 

DIST=GSDIS (LONASSET (ASSET) , LATASSET ( ASSET ) , 
Sc LON(Q) , LAT(Q) ) 

END IF 



* COMPUTE PROBABILITIES (INDEPENDENCE ASSUMEDO 
IF ( DIST .LE. 

Sc RANGE (ASSETTY PE (ASSET) , ROUTEMETHOD ( R ) ) ) THEN 

A = 1 . 0 - I PD ( ASSETTYPE ( ASSET ) , ROUTEMETHOD ( R) ) 



Sc 

Sc 

Sc 



420 



Sc 



ELSE 

IF (RANGE (ASSETTYPE (ASSET) , ROUTEMETHOD ( R ) ) . EQ . 0.0) 

RANGE (ASSETTYPE (ASSET) , ROUTEMETHOD (R) ) = 0.0001 
FCTR= (RANGE (ASSETTYPE (ASSET) , ROUTEMETHOD ( R) ) -DIST) / 
RANGE (ASSETTYPE (ASSET) , ROUTEMETHOD { R ) ) 

A = 1.0 - I PD { ASSETTYPE ( ASSET ), ROUTEMETHOD (R) ) 

* EXP(FCTR) 

END IF 

TEMPPI =1.0- (A * (1.0- TEMPPI ) ) 

CONTINUE 

DO 430 N=1 , NUMPHASES 

PROB INTERDICT ( R, N) = 1.0 -{1.0 - TEMPPI ) * 

(1.0 - PROBINTERDICT ( R , N ) ) 



430 CONTINUE 

P = Q 

Q = NEXT { R, P) 

IF { Q .NE. 0) GO TO 410 
440 CONTINUE 
RETURN 
END 






REAL FUNCTION GSDIS( XI, Yl, X2 , Y2 ) 
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FUNCTION COMPUTES THE GREAT CIRCLE DISTANCE IN MILES 
BETWEEN TWO POINTS EXPRESSED IN LATITUDE AND LONGITUDE 



REAL *4 
REAL *4 



XI, Y1,X2,Y2 

RARC , RLAT1 , RLAT2 , RLON1 , RLON2 , RPD 



--LOCAL DATA 
DATA 



INITIALIZATION 

RPD /0.1745329433017307E-01/ 



--GREAT CIRCLE DISTANCE 

INPUTS ARE TWO POINTS EXPRESSED AS {W. LONGITUDE, N. LATITUDE) 
-- CONVERT DEGREES TO RADIANS 
RLON1 = RPD * XI 

RLAT1 = RPD * Y1 

RLON2 = RPD * X2 

RLAT2 = RPD * Y2 

--SPHERICAL ARC LENGTH IN RADIANS 

ARGUMENT = { SIN { RLAT1 ) *SIN ( RLAT2 ) 

& +COS { RLAT1 ) *COS { RLAT2 ) *COS { RLON1 -RLON2 ) ) 

IF (ARGUMENT .GT. 1.0) ARGUMENT =1.0 
IF (ARGUMENT . LT . -1.0) ARGUMENT = -1.0 
RARC = ACOS( ARGUMENT) 

--ARC LENGTH IN STATUTE MILES WITH ROUGH CORRECTION FOR OBLATENESS 
OBLATENESS CONTRIBUTES LESS THAN 0.25% CORRECTION. 

IN U.S., DARC IS APPROXIMATELY 3959.U*RARC 

GSDIS = 3959 . 0E0*RARC 



***** RANDOM NUMBER GENERATING SUBROUTINES ******************** 
ADAPTED FROM LEWIS , ORAV, AND URIBE ENHANCED SIMULATION AND 
STATISTICS PACKAGE [Ref. 27] 

THIS PROGRAM WILL GENERATE A VECTOR OF NORMAL RANDOM VARIABLES 
ACCORDING TO THE SINE-COSINE METHOD 

SUBROUTINE LNORPC ( ISEED, A, N) 

INTEGER N, I , IND, ISEED 
REAL A ( N ) , U ( 2 ) ,S,W,XSTAR 
DOUBLE PRECISION PI 
DATA PI/3 .14159265358979 DO/ 

IND = 1 
DO 100 1=1, N 

IND = -IND 

IF (IND.GE.O) GOTO 20 
10 CALL LRNDPC ( ISEED, U , 2 , 0 , 0 ) 

S = SQRT ( - 2 * ALOG ( U ( 1 ) ) ) 

W = 2 * PI *U ( 2 ) 

XSTAR = S*COS(W) 

A ( I ) = S*SIN (W) 

GOTO 100 

20 A ( I ) = XSTAR 

100 CONTINUE 
RETURN 
END 



RETURN 

END 



SUBROUTINE LRNDPC 


( ISEED, U 


, N) 


INTEGER 


N, 1, ISEED 


REAL 


U ( N ) 




DOUBLE PRECISION 
D3 1M1 =2**31 - 1 
D31 =2**31 


D31M1, 


DSEED 


DATA D31M1/2147483647 . DO/ 
DATA D31 /2147483648 .DO/ 
DSEED = DBLE( ISEED) 




DO 5 1 = 1, N 







C DSEED = DMOD (950706376 . DO * DSEED, D3 1M1 ) 

DSEED = DMOD ( 16807 . DO *DSEED, D3 1M1 ) 

U(I) = DSEED / D31 
5 CONTINUE 

ISEED = I NT (DSEED) 

RETURN 

END 



★ ★★★★★* 



COMMON . FOR 



C SET PAMETERS : 
INTEGER* 4 

Sc 

Sc 

L 



LONG PAST, RECENTPAST, MAXDAYS , 

MAXMETHODS, MAXPHASES, MAXROUTES , MAXSH I PMENTS , 
MAXASSETS , MAXASSETTYPES , MAXPOINTS , 

MAXPTS 



PARAMETER 

Sc 

Sc 

Sc 

Sc 



( LONG P AST = 12 0 , RECENT PAST = 2 0 , 

MAXDAYS=73 0 , 

MAXMETHODS =3 , MAXPHASES=12 , MAX ROUTES =100 , 
MAXSH I PMENTS = 6000, MAXASSETS = 1 0 0 , 
MAXASSETTYPES =4 , MAXPOINTS=50 , MAXPTS=8 } 



C COMMON VARIABLES: 



C 



REAL 

REAL 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 

Sc 



VARIABLES : 

AMOUNT ATTEMPTED { MAXROUTES) , 
AMOUNTSUCCEEDED ( MAXROUTES) , 

ATTEMPTS (MAXROUTES) , 

ATTEMPTSBYPHASE (MAXROUTES , MAXPHASES ) , 
CAPCOST (MAXMETHODS) , 

CAPACITY (MAXMETHODS) , 

DAILYAMOUNT, 

DRUGCOST , EX PT ABLE ( 0 : LONGPAST ) , 

EXSHI PMENT INTRVL , 

EXSH I PMENTS I ZE , 

FAI LURECOSTS , KMEMORY , MEMORYVALUE , 
NEXTEVENT , NEXTSH I PMENT , 

PROBINTERDICT (MAXROUTES , MAXMETHODS ) , 

R I SKCOMP (MAXMETHODS) , 

RISKCOMPEXP (MAXMETHODS) , 

ROUTECOST (MAXROUTES) , 

SUCCESSCOSTS, 

SUCCESSES (MAXROUTES) , 
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tr> & 



Zc SUCCESSESBYPHASE (MAXROUTES , MAXPHASES ) , 

Z LAT( MAX POINTS) ( LON (MAX POINTS ) , 

LAT ASSET (MAXASSETS) , LONASSET (MAXASSETS) , 
DELR , DELRMIN, DELRLOW, 

BETA, ALPHA, RANGE (MAXASSETTYPES , MAXMETHODS ) , 
Zc I PD ( MAXASSETT Y PES , MAXMETHODS ), GUESS 



C 



Zc 

Zc 

Zc 

Zc 

Zc 

& 

Zc 

Zc 

Zc 

& 

Zc 

Zc 

Zc 



INTEGER VARIABLES: 

INTEGER* 4 CURRENTPHASE { 0 : MAXDAYS ) , DAYNOW, ENDTIME, NTRIAL, 

NEXTEVENTTYPE, NUMMETHODS , NUMPHASES, 

NUMROUTES, NUMSHIPMENT, NUMTRIALS , 

PASTSH I PMENTS { MAXROUTES , 0 : LONG PAST ) , 
PASTFAILURES (MAXROUTES, 0:LONGPAST) , 

ROUTEMETHOD (MAXROUTES) ,RUNIN, 

THISPHASE, 

PEACESH I PMENTS (MAXROUTES, 0:LONGPAST) , 
PEACEFAILURES (MAXROUTES, OiLONGPAST) , 

NUMASSETS, ROUTELOC (MAXROUTES) , 

ASSETLOC (MAXASSETS) , 

ASSETTY PE ( MAXASSETS ) , NUMASSETTYPES , 

S, NEXT (MAXROUTES, MAXPOINTS) , 

FIRST (MAXROUTES) , BIGN 



C INTEGER VARIABLES: 

INTEGER *4 IX, ISEED 



C CHARACTER VARIABLES: 

CHARACTER* 10 METHODNAME (MAXMETHODS ) 



CHARACTER * 1 2 ROUTENAME ( MAXROUTES ) 

CHARACTER * 1 5 ASSETNAME ( MAX ASSETS ) 



COMMON 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 

Zc 



AMOUNTATTEMPTED, AMOUNTSUCCEEDED, ATTEMPTS, 
ATTEMPTSBYPHASE, 

CAPCOST , 

CAPACITY, 

DAILY AMOUNT, 

DRUGCOST, EXPTABLE, 

EXSH I PMENT INTRVL, 

EXSH I PMENTS I ZE, 

FAILURECOSTS , KMEMORY, MEMORYVALUE, 
NEXTEVENT , NEXTSHI PMENT , 
PROBINTERDICT, 

RISKCOMP, 

RISKCOMPEXP , 

ROUTECOST, 

SUCCESSCOSTS, 

SUCCESSES, 

SUCCESSESBYPHASE, 

LAT, LON, 

LATASSET, LONASSET, 

DELR, DELRMIN, DELR LOW, 

BETA, ALPHA, RANGE, IPD, GUESS, 
CURRENTPHASE, DAYNOW, ENDTIME, NTRIAL, 
NEXTEVENTTYPE, NUMMETHODS , NUMPHASES , 
NUMROUTES , NUMSHI PMENT , NUMTR I ALS , 
PASTSH I PMENTS, 

PASTFAILURES, 

ROUTEMETHOD, RUNIN, 

THISPHASE, 
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c^> c^> 



& PEACESH I PMENTS # 

Sc PEACEFAILURES , 

Sc NUMASSETS , ROUTELOC , 

Sc ASSETLOC , 

& ASSETTYPE, NUMASSETTYPES, 

& S , NEXT , FIRST, BIGN, 

IX, ISEED, 

METHODNAME, 

ROUTENAME, ASSETNAME 
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